1 Setup

1.1 Load libraries

library(survival)
library(piecewiseSEM)
source("spat1f/init_analysis.R")
Data last updated on 2024-05-03

1.2 Handy functions not in spat1f

trend_by_pre_wt <- function(model, term){
  emmeans::emtrends(model, var = "cat_pre_wt_log_scale", specs = term) %>% 
  pairs()
}

contrast_by_pre_wt <- function(model, term, type = c("pairwise", "trt.vs.ctrl"), at = c(-2, 2)){
  type <- match.arg(type)
  f <- as.formula(sprintf("%s ~ %s | cat_pre_wt_log_scale", type, term))
  emmeans::emmeans(model, 
                 f, 
                 var = term, 
                 at = list(cat_pre_wt_log_scale = at))
}

check_mod <- function(model){
 DHARMa::simulateResiduals(model) %>% 
    plot()
  performance::posterior_predictive_check(model)
}

# Get posterior draws
epred_draws <- function(model, newdata){
  res <- brms::posterior_epred(model, 
                  newdata = newdata, 
                  allow_new_levels = TRUE, 
                  re_formula = NA) %>% 
    t()
  colnames(res) <- paste0("draw_", seq_len(ncol(res)))
  cbind(newdata,res) %>% 
    gather(key = draw, value = val, contains("draw_"))
}

1.3 Quick cleaning

# Main data.frame for analyses
d <- ref_data %>% 
  filter(error == 0) %>% 
  mutate(
    cat_pre_wt_log = log(cat_pre_wt),
    cat_pre_wt_log_scale = as.numeric(scale(log(cat_pre_wt))),
    cat_size = ifelse(cat_pre_wt < median(cat_pre_wt), "small", "big"),
    has_var = ifelse(var_trt != "constant", 1, 0),
    mean_toxic_conc_scale = as.numeric(scale(mean_toxic_conc)),
    mean_toxic_scale = as.numeric(scale(mean_toxic)),
    area_herb_log_scale = as.numeric(scale(log(area_herb+1))),
    var_toxic_12_scale = as.numeric(scale(var_toxic_12)),
    on_toxic_logit_scale = as.numeric(scale(adjust_prop(on_toxic, 
                                                        trans = "emp", 
                                                        nudge.method = "none", 
                                                        na.action = "ignore"))),
    sl_mean_obs_log_scale = as.numeric(scale(log(sl_mean_obs))),
    cat_pre_wt_log_scale_sq = as.numeric(scale(log(cat_pre_wt)^2)),
    RGR_scale = as.numeric(scale(RGR)), 
    sl_mean_pred1 = shape_1 * scale_1, # Exploration state & on more toxic diet 
    sl_mean_pred2 = shape_2 * scale_2, # Resting/feeding state & on more toxic diet
    sl_mean_pred3 = shape_3 * scale_3, # Exploration state & on less toxic diet
    sl_mean_pred4 = shape_4 * scale_4, # Resting/feeding state & on less toxic diet
    k1 = scale_1 / scale_3, # state 1 (exploration)
    k2 = scale_2 / scale_4, # state 2 (resting/feeding)
    prop_explore_logit = qlogis(prop_explore),
    prop_explore_logit_scale = as.numeric(scale(prop_explore_logit))
  )

# Subset of data.frame for SEM
d2 <- d %>% 
  filter(
    var_trt != "constant"
  ) %>% 
  filter(
    !is.na(area_herb_log_scale) & 
      !is.na(var_toxic_12_scale) & 
      !is.na(mean_toxic_conc_scale) &
      !is.na(sl_mean_obs) &
      !is.na(prop_explore) & 
      !is.na(RGR)) %>% 
  mutate(
    var_high = ifelse(var_trt == "high_var", 1, 0),
    beta_red = ifelse(as.character(beta) == 5, 1, 0),
    beta_white = ifelse(as.character(beta) == 0, 1, 0),
    beta_red_cat_pre_wt_log_scale = beta_red * cat_pre_wt_log_scale,
    beta_white_cat_pre_wt_log_scale = beta_white * cat_pre_wt_log_scale,
    var_high_cat_pre_wt_log_scale = var_high * cat_pre_wt_log_scale,
    beta_numeric_scale = as.numeric(scale(as.numeric(beta)))
  )

# Subset of data.frame for more detailed analyses
d3 <- d %>% 
  filter(
    var_trt != "constant"
  ) %>% 
  filter_at(vars(contains("shape_"), contains("scale_[0-9]"), contains("sl_mean_pred")), 
            all_vars(. > 0)) %>% 
  mutate(
    var_high = ifelse(var_trt == "high_var", 1, 0),
    beta_red = ifelse(as.character(beta) == 5, 1, 0),
    beta_white = ifelse(as.character(beta) == 0, 1, 0),
    beta_red_cat_pre_wt_log_scale = beta_red * cat_pre_wt_log_scale,
    beta_white_cat_pre_wt_log_scale = beta_white * cat_pre_wt_log_scale,
    var_high_cat_pre_wt_log_scale = var_high * cat_pre_wt_log_scale,
    beta_numeric_scale = as.numeric(scale(as.numeric(beta)))
  )

2 Herbivore Performance

2.1 RGR

2.1.1 Model summary

rgr_m <- glmmTMB(
  RGR ~ 
    (var_trt + beta) * cat_pre_wt_log_scale +
    I(cat_pre_wt_log_scale^2) +  # residual plot show significant non-linearity
    (1|session_id),
  data = d %>% 
    filter(
      var_trt != "constant"
    ) 
); summary(rgr_m)
 Family: gaussian  ( identity )
Formula:          
RGR ~ (var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +  
    (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
  -879.2   -849.4    450.6   -901.2      100 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev.
 session_id (Intercept) 3.979e-06 0.001995
 Residual               1.603e-05 0.004004
Number of obs: 111, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 1.6e-05 

Conditional model:
                                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          0.0162504  0.0013457  12.075  < 2e-16 ***
var_trtlow_var                      -0.0034872  0.0007740  -4.506 6.62e-06 ***
beta0                               -0.0003065  0.0009664  -0.317  0.75112    
beta5                               -0.0013545  0.0009173  -1.477  0.13978    
cat_pre_wt_log_scale                -0.0048759  0.0011991  -4.066 4.78e-05 ***
I(cat_pre_wt_log_scale^2)           -0.0021953  0.0006731  -3.262  0.00111 ** 
var_trtlow_var:cat_pre_wt_log_scale  0.0025873  0.0007897   3.276  0.00105 ** 
beta0:cat_pre_wt_log_scale           0.0011924  0.0009491   1.256  0.20897    
beta5:cat_pre_wt_log_scale           0.0024731  0.0009710   2.547  0.01087 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(rgr_m, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: RGR
                                Chisq Df Pr(>Chisq)    
(Intercept)                  145.8141  1  < 2.2e-16 ***
var_trt                       20.3008  1  6.617e-06 ***
beta                           2.3963  2   0.301758    
cat_pre_wt_log_scale          16.5348  1  4.777e-05 ***
I(cat_pre_wt_log_scale^2)     10.6374  1   0.001108 ** 
var_trt:cat_pre_wt_log_scale  10.7334  1   0.001052 ** 
beta:cat_pre_wt_log_scale      6.4936  2   0.038899 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.1.2 Post-hoc contrasts

trend_by_pre_wt(rgr_m, "var_trt")
 contrast           estimate      SE  df t.ratio p.value
 high_var - low_var -0.00259 0.00079 100  -3.276  0.0014

Results are averaged over the levels of: beta 
contrast_by_pre_wt(rgr_m, "var_trt")
$emmeans
cat_pre_wt_log_scale = -2:
 var_trt     emmean      SE  df  lower.CL upper.CL
 high_var  0.014224 0.00355 100  0.007178  0.02127
 low_var   0.005562 0.00312 100 -0.000635  0.01176

cat_pre_wt_log_scale =  2:
 var_trt     emmean      SE  df  lower.CL upper.CL
 high_var -0.000393 0.00208 100 -0.004512  0.00373
 low_var   0.001295 0.00238 100 -0.003421  0.00601

Results are averaged over the levels of: beta 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast           estimate      SE  df t.ratio p.value
 high_var - low_var  0.00866 0.00181 100   4.790  <.0001

cat_pre_wt_log_scale =  2:
 contrast           estimate      SE  df t.ratio p.value
 high_var - low_var -0.00169 0.00171 100  -0.988  0.3256

Results are averaged over the levels of: beta 
trend_by_pre_wt(rgr_m, "beta")
 contrast         estimate       SE  df t.ratio p.value
 (beta-5) - beta0 -0.00119 0.000949 100  -1.256  0.4232
 (beta-5) - beta5 -0.00247 0.000971 100  -2.547  0.0329
 beta0 - beta5    -0.00128 0.000940 100  -1.363  0.3644

Results are averaged over the levels of: var_trt 
P value adjustment: tukey method for comparing a family of 3 estimates 
contrast_by_pre_wt(rgr_m, "beta")
$emmeans
cat_pre_wt_log_scale = -2:
 beta    emmean      SE  df  lower.CL upper.CL
 -5    0.012890 0.00382 100  0.005313  0.02047
 0     0.010199 0.00333 100  0.003591  0.01681
 5     0.006589 0.00319 100  0.000255  0.01292

cat_pre_wt_log_scale =  2:
 beta    emmean      SE  df  lower.CL upper.CL
 -5   -0.001439 0.00226 100 -0.005922  0.00304
 0     0.000639 0.00231 100 -0.003940  0.00522
 5     0.002153 0.00257 100 -0.002946  0.00725

Results are averaged over the levels of: var_trt 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast         estimate      SE  df t.ratio p.value
 (beta-5) - beta0  0.00269 0.00221 100   1.218  0.4456
 (beta-5) - beta5  0.00630 0.00223 100   2.820  0.0158
 beta0 - beta5     0.00361 0.00211 100   1.711  0.2062

cat_pre_wt_log_scale =  2:
 contrast         estimate      SE  df t.ratio p.value
 (beta-5) - beta0 -0.00208 0.00205 100  -1.015  0.5689
 (beta-5) - beta5 -0.00359 0.00206 100  -1.746  0.1935
 beta0 - beta5    -0.00151 0.00211 100  -0.719  0.7529

Results are averaged over the levels of: var_trt 
P value adjustment: tukey method for comparing a family of 3 estimates 

2.1.3 Misc stats

r2(rgr_m, tolerance = 10^-10)
# R2 for Mixed Models

  Conditional R2: 0.555
     Marginal R2: 0.444

2.1.4 Check model

check_mod(rgr_m)

2.1.5 Plot

g1 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","beta")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) + 
  geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) + 
  geom_line(aes(color = beta), linewidth = 2) + 
  geom_point(
    data = rgr_m$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = beta),
    size = 3, alpha = 0.8
  ) + 
  theme_bw(base_size = 15) + 
  theme(legend.position = "top") +
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  scale_color_brewer(type = "qual", aesthetics = c("fill","color")) + 
  labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
       color = expression(beta), fill = expression(beta))

g2 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","var_trt")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) + 
  geom_ribbon(aes(ymax = upper, ymin = lower, fill = var_trt), alpha = 0.2) + 
  geom_line(aes(color = var_trt), linewidth = 2) + 
  geom_point(
    data = rgr_m$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = var_trt),
    size = 3, alpha = 0.8
  ) + 
  theme_bw(base_size = 15) + 
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  theme(legend.position = "top") +
  scale_color_discrete(type = c("navy","skyblue"), 
                       aesthetics = c("fill","color"), 
                       label = c("High", "Low")) + 
  labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
       color = "Variation", fill = "Variation")


ggarrange(g1, g2 + labs(y = ""))

2.2 Time to pupation

2.2.1 Model summary

pupt_m_full <- glmmTMB(
  time_to_pupation ~ 
    (var_trt + beta) * cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = Gamma(link = "log"),
); summary(pupt_m_full)
 Family: Gamma  ( log )
Formula:          time_to_pupation ~ (var_trt + beta) * cat_pre_wt_log_scale +  
    I(cat_pre_wt_log_scale^2) + (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   338.1    366.3   -158.0    316.1       85 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 2.839e-11 5.328e-06
Number of obs: 96, groups:  session_id, 5

Dispersion estimate for Gamma family (sigma^2): 0.0254 

Conditional model:
                                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          2.180302   0.036853   59.16  < 2e-16 ***
var_trtlow_var                       0.102145   0.033998    3.00  0.00266 ** 
beta0                               -0.001845   0.040823   -0.05  0.96395    
beta5                                0.059543   0.041181    1.45  0.14821    
cat_pre_wt_log_scale                -0.428662   0.032643  -13.13  < 2e-16 ***
I(cat_pre_wt_log_scale^2)           -0.081267   0.019069   -4.26 2.03e-05 ***
var_trtlow_var:cat_pre_wt_log_scale  0.026518   0.038876    0.68  0.49516    
beta0:cat_pre_wt_log_scale          -0.030363   0.043195   -0.70  0.48210    
beta5:cat_pre_wt_log_scale          -0.126801   0.047095   -2.69  0.00709 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(pupt_m_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: time_to_pupation
                                 Chisq Df Pr(>Chisq)    
(Intercept)                  3500.2207  1  < 2.2e-16 ***
var_trt                         9.0264  1   0.002661 ** 
beta                            2.7351  2   0.254736    
cat_pre_wt_log_scale          172.4448  1  < 2.2e-16 ***
I(cat_pre_wt_log_scale^2)      18.1614  1  2.029e-05 ***
var_trt:cat_pre_wt_log_scale    0.4653  1   0.495162    
beta:cat_pre_wt_log_scale       7.5589  2   0.022836 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pupt_m <- glmmTMB(
  time_to_pupation ~ 
    var_trt + beta * cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = Gamma(link = "log"),
); summary(pupt_m)
 Family: Gamma  ( log )
Formula:          
time_to_pupation ~ var_trt + beta * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +  
    (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   336.5    362.2   -158.3    316.5       86 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 2.159e-11 4.647e-06
Number of obs: 96, groups:  session_id, 5

Dispersion estimate for Gamma family (sigma^2): 0.0255 

Conditional model:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)                 2.180555   0.036937   59.03  < 2e-16 ***
var_trtlow_var              0.107447   0.033197    3.24  0.00121 ** 
beta0                      -0.003578   0.040842   -0.09  0.93019    
beta5                       0.059536   0.041271    1.44  0.14915    
cat_pre_wt_log_scale       -0.420147   0.030221  -13.90  < 2e-16 ***
I(cat_pre_wt_log_scale^2)  -0.083365   0.018881   -4.42 1.01e-05 ***
beta0:cat_pre_wt_log_scale -0.026011   0.042849   -0.61  0.54383    
beta5:cat_pre_wt_log_scale -0.126834   0.047183   -2.69  0.00719 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(pupt_m, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: time_to_pupation
                              Chisq Df Pr(>Chisq)    
(Intercept)               3485.0791  1  < 2.2e-16 ***
var_trt                     10.4759  1   0.001209 ** 
beta                         2.8053  2   0.245947    
cat_pre_wt_log_scale       193.2743  1  < 2.2e-16 ***
I(cat_pre_wt_log_scale^2)   19.4946  1  1.009e-05 ***
beta:cat_pre_wt_log_scale    7.6890  2   0.021397 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.2.2 Post-hoc contrasts

contrast_by_pre_wt(pupt_m, "var_trt")
$emmeans
cat_pre_wt_log_scale = -2:
 var_trt  emmean     SE  df asymp.LCL asymp.UCL
 high_var  2.808 0.0769 Inf     2.657      2.96
 low_var   2.915 0.0818 Inf     2.755      3.08

cat_pre_wt_log_scale =  2:
 var_trt  emmean     SE  df asymp.LCL asymp.UCL
 high_var  0.924 0.0672 Inf     0.792      1.06
 low_var   1.031 0.0725 Inf     0.889      1.17

Results are averaged over the levels of: beta 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast           estimate     SE  df z.ratio p.value
 high_var - low_var   -0.107 0.0332 Inf  -3.237  0.0012

cat_pre_wt_log_scale =  2:
 contrast           estimate     SE  df z.ratio p.value
 high_var - low_var   -0.107 0.0332 Inf  -3.237  0.0012

Results are averaged over the levels of: beta 
Results are given on the log (not the response) scale. 
trend_by_pre_wt(pupt_m, "beta")
 contrast         estimate     SE  df z.ratio p.value
 (beta-5) - beta0    0.026 0.0428 Inf   0.607  0.8163
 (beta-5) - beta5    0.127 0.0472 Inf   2.688  0.0197
 beta0 - beta5       0.101 0.0473 Inf   2.132  0.0834

Results are averaged over the levels of: var_trt 
P value adjustment: tukey method for comparing a family of 3 estimates 
contrast_by_pre_wt(pupt_m, "beta")
$emmeans
cat_pre_wt_log_scale = -2:
 beta emmean     SE  df asymp.LCL asymp.UCL
 -5    2.741 0.0941 Inf     2.557      2.93
 0     2.790 0.0971 Inf     2.599      2.98
 5     3.054 0.1088 Inf     2.841      3.27

cat_pre_wt_log_scale =  2:
 beta emmean     SE  df asymp.LCL asymp.UCL
 -5    1.061 0.0845 Inf     0.895      1.23
 0     1.005 0.0805 Inf     0.847      1.16
 5     0.866 0.0917 Inf     0.687      1.05

Results are averaged over the levels of: var_trt 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast         estimate     SE  df z.ratio p.value
 (beta-5) - beta0  -0.0484 0.1018 Inf  -0.476  0.8828
 (beta-5) - beta5  -0.3132 0.1133 Inf  -2.763  0.0158
 beta0 - beta5     -0.2648 0.1136 Inf  -2.331  0.0516

cat_pre_wt_log_scale =  2:
 contrast         estimate     SE  df z.ratio p.value
 (beta-5) - beta0   0.0556 0.0875 Inf   0.636  0.8005
 (beta-5) - beta5   0.1941 0.0915 Inf   2.122  0.0854
 beta0 - beta5      0.1385 0.0928 Inf   1.492  0.2945

Results are averaged over the levels of: var_trt 
Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 

2.2.3 Misc stats

r2(pupt_m, tolerance = 10^-10)
Random effect variances not available. Returned R2 does not account for random effects.
# R2 for Mixed Models

  Conditional R2: NA
     Marginal R2: 0.881

2.2.4 Check model

check_mod(pupt_m)

2.2.5 Plot

g3 <- marginal_effects(pupt_m, terms = c("cat_pre_wt_log_scale","beta")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) + 
  geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) + 
  geom_line(aes(color = beta), linewidth = 2) + 
  geom_point(
    data = pupt_m$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = time_to_pupation, color = beta),
    size = 3, alpha = 0.8
  ) + 
  theme_bw(base_size = 15) + 
  theme(legend.position = "top") +
  scale_y_continuous(trans = "log2") +
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  scale_color_brewer(type = "qual", aesthetics = c("fill","color")) + 
  labs(x = "Cat pre-weight (g)", y = "Time to pupation (days)",
       color = expression(beta), fill = expression(beta))

g4 <- marginal_effects(pupt_m, terms = c("var_trt")) %>% 
  ggplot(aes(x = var_trt, y = yhat)) + 
  geom_point(
    data = pupt_m$frame,
    aes(x = var_trt, y = time_to_pupation, color = var_trt),
    position = position_jitter(width = 0.2),
    size = 3, alpha = 0.8
  ) + 
  geom_pointrange(
    aes(ymax = upper, ymin = lower), 
    color = "black",
    size = 1, linewidth = 2, shape = 1,
  ) +
  theme_bw(base_size = 15) + 
  scale_y_continuous(trans = "log2") +
  theme(legend.position = "top") +
  scale_color_discrete(type = c("navy","skyblue"), 
                       aesthetics = c("color"), 
                       label = c("High", "Low")) + 
  labs(x = "Variation", y = "Time to pupation (days)",
       color = "Variation")


ggarrange(g3, g4 + labs(y = ""))

2.3 Eclosure

2.3.1 Model summary

eclosure_m_full <- glmmTMB(
  eclosed ~ 
    (var_trt + beta) * cat_pre_wt_log_scale + 
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = binomial(link = "logit"),
); summary(eclosure_m_full)
 Family: binomial  ( logit )
Formula:          
eclosed ~ (var_trt + beta) * cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   185.3    210.9    -83.6    167.3      119 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.1661   0.4076  
Number of obs: 128, groups:  session_id, 5

Conditional model:
                                    Estimate Std. Error z value Pr(>|z|)
(Intercept)                          0.53530    0.41050   1.304    0.192
var_trtlow_var                       0.01707    0.37507   0.046    0.964
beta0                               -0.21544    0.46265  -0.466    0.641
beta5                               -0.50749    0.45596  -1.113    0.266
cat_pre_wt_log_scale                 0.10452    0.36089   0.290    0.772
var_trtlow_var:cat_pre_wt_log_scale  0.17526    0.36489   0.480    0.631
beta0:cat_pre_wt_log_scale           0.24778    0.43553   0.569    0.569
beta5:cat_pre_wt_log_scale           0.20159    0.44699   0.451    0.652
car::Anova(eclosure_m_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: eclosed
                              Chisq Df Pr(>Chisq)
(Intercept)                  1.7004  1     0.1922
var_trt                      0.0021  1     0.9637
beta                         1.2502  2     0.5352
cat_pre_wt_log_scale         0.0839  1     0.7721
var_trt:cat_pre_wt_log_scale 0.2307  1     0.6310
beta:cat_pre_wt_log_scale    0.3659  2     0.8328
eclosure_m <- glmmTMB(
  eclosed ~ 
    (var_trt + beta) + cat_pre_wt_log_scale + 
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = binomial(link = "logit"),
); summary(eclosure_m)
 Family: binomial  ( logit )
Formula:          
eclosed ~ (var_trt + beta) + cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   179.9    197.0    -84.0    167.9      122 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.1676   0.4094  
Number of obs: 128, groups:  session_id, 5

Conditional model:
                      Estimate Std. Error z value Pr(>|z|)
(Intercept)           0.551286   0.413578   1.333    0.183
var_trtlow_var        0.005064   0.373402   0.014    0.989
beta0                -0.237692   0.460689  -0.516    0.606
beta5                -0.523521   0.456587  -1.147    0.252
cat_pre_wt_log_scale  0.330169   0.220049   1.500    0.134
car::Anova(eclosure_m, type = "II")
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: eclosed
                      Chisq Df Pr(>Chisq)
var_trt              0.0002  1     0.9892
beta                 1.3212  2     0.5165
cat_pre_wt_log_scale 2.2513  1     0.1335

2.3.2 Post-hoc contrasts

2.3.3 Misc stats

r2(eclosure_m, tolerance = 10^-10)
# R2 for Mixed Models

  Conditional R2: 0.092
     Marginal R2: 0.046

2.3.4 Check model

check_mod(eclosure_m)

2.3.5 Plot

2.4 Survival

2.4.1 Model summary

surv_m_full <- coxph(
  Surv(surv_time, observed_dead) ~ 
    (var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
    strata(session_id),
  data = d %>% 
    filter(var_trt != "constant") %>% 
    mutate(
      beta = as.factor(as.character(beta))
    ),
); summary(surv_m_full)
Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) * 
    cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id), 
    data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))

  n= 128, number of events= 105 

                                       coef exp(coef) se(coef)      z Pr(>|z|)
var_trtlow_var                      -0.1115    0.8945   0.2124 -0.525 0.599727
beta0                               -0.0896    0.9143   0.2579 -0.347 0.728306
beta5                               -0.2521    0.7772   0.2605 -0.968 0.333144
cat_pre_wt_log_scale                 0.9641    2.6224   0.2556  3.772 0.000162
I(cat_pre_wt_log_scale^2)            0.2450    1.2776   0.1447  1.693 0.090408
var_trtlow_var:cat_pre_wt_log_scale -0.1626    0.8499   0.1998 -0.814 0.415661
beta0:cat_pre_wt_log_scale          -0.2230    0.8001   0.2309 -0.966 0.334053
beta5:cat_pre_wt_log_scale          -0.4789    0.6195   0.2584 -1.854 0.063807
                                       
var_trtlow_var                         
beta0                                  
beta5                                  
cat_pre_wt_log_scale                ***
I(cat_pre_wt_log_scale^2)           .  
var_trtlow_var:cat_pre_wt_log_scale    
beta0:cat_pre_wt_log_scale             
beta5:cat_pre_wt_log_scale          .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                                    exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var                         0.8945     1.1179    0.5899     1.356
beta0                                  0.9143     1.0937    0.5515     1.516
beta5                                  0.7772     1.2867    0.4665     1.295
cat_pre_wt_log_scale                   2.6224     0.3813    1.5892     4.328
I(cat_pre_wt_log_scale^2)              1.2776     0.7827    0.9621     1.696
var_trtlow_var:cat_pre_wt_log_scale    0.8499     1.1766    0.5746     1.257
beta0:cat_pre_wt_log_scale             0.8001     1.2499    0.5089     1.258
beta5:cat_pre_wt_log_scale             0.6195     1.6142    0.3734     1.028

Concordance= 0.637  (se = 0.047 )
Likelihood ratio test= 22  on 8 df,   p=0.005
Wald test            = 22.31  on 8 df,   p=0.004
Score (logrank) test = 24.57  on 8 df,   p=0.002
drop1(surv_m_full, test = "Chisq")
Single term deletions

Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) * cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + strata(session_id)
                             Df    AIC    LRT Pr(>Chi)  
<none>                          459.50                  
I(cat_pre_wt_log_scale^2)     1 460.40   2.90  0.08841 .
strata(session_id)            0 787.87 328.37           
var_trt:cat_pre_wt_log_scale  1 458.16   0.67  0.41395  
beta:cat_pre_wt_log_scale     2 458.98   3.48  0.17532  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
surv_m <- coxph(
  Surv(surv_time, observed_dead) ~ 
    (var_trt + beta) + cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + 
    strata(session_id),
  data = d %>% 
    filter(var_trt != "constant") %>% 
    mutate(
      beta = as.factor(as.character(beta))
    ),
); summary(surv_m)
Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) + 
    cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id), 
    data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))

  n= 128, number of events= 105 

                              coef exp(coef) se(coef)      z Pr(>|z|)    
var_trtlow_var            -0.07077   0.93167  0.20845 -0.340 0.734212    
beta0                     -0.07410   0.92858  0.25808 -0.287 0.774026    
beta5                     -0.17722   0.83759  0.25804 -0.687 0.492206    
cat_pre_wt_log_scale       0.65947   1.93376  0.18164  3.631 0.000283 ***
I(cat_pre_wt_log_scale^2)  0.29584   1.34425  0.14124  2.095 0.036205 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                          exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var               0.9317     1.0733    0.6192     1.402
beta0                        0.9286     1.0769    0.5599     1.540
beta5                        0.8376     1.1939    0.5051     1.389
cat_pre_wt_log_scale         1.9338     0.5171    1.3545     2.761
I(cat_pre_wt_log_scale^2)    1.3443     0.7439    1.0192     1.773

Concordance= 0.602  (se = 0.047 )
Likelihood ratio test= 18.08  on 5 df,   p=0.003
Wald test            = 18.29  on 5 df,   p=0.003
Score (logrank) test = 19.42  on 5 df,   p=0.002
drop1(surv_m, test = "Chisq")
Single term deletions

Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) + cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + strata(session_id)
                          Df    AIC    LRT  Pr(>Chi)    
<none>                       457.42                     
var_trt                    1 455.53   0.12 0.7341150    
beta                       2 453.89   0.48 0.7883535    
cat_pre_wt_log_scale       1 469.05  13.63 0.0002226 ***
I(cat_pre_wt_log_scale^2)  1 459.84   4.42 0.0355376 *  
strata(session_id)         0 784.70 327.28              
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.4.2 Post-hoc contrasts

2.4.3 Misc stats

r2(surv_m, tolerance = 10^-10)
  Nagelkerke's R2: 0.135

2.4.4 Check model

cox.zph(surv_m)
                          chisq df    p
var_trt                    1.60  1 0.21
beta                       1.43  2 0.49
cat_pre_wt_log_scale       2.61  1 0.11
I(cat_pre_wt_log_scale^2)  1.05  1 0.31
GLOBAL                     5.53  5 0.35

2.4.5 Plot

3 Structural Equation Model

3.1 Sub-models

3.1.1 RGR

subm_rgr <- glmmTMB(
  RGR_scale ~ 
    cat_pre_wt_log_scale + 
    cat_pre_wt_log_scale_sq + 
    sl_mean_obs_log_scale +
    prop_explore_logit_scale *  
    var_toxic_12_scale + 
    mean_toxic_conc_scale + 
    area_herb_log_scale + 
    beta_numeric_scale + 
    (1|session_id),
  data = d2
); summary(subm_rgr)
 Family: gaussian  ( identity )
Formula:          RGR_scale ~ cat_pre_wt_log_scale + cat_pre_wt_log_scale_sq +  
    sl_mean_obs_log_scale + prop_explore_logit_scale * var_toxic_12_scale +  
    mean_toxic_conc_scale + area_herb_log_scale + beta_numeric_scale +  
    (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   119.5    149.0    -47.7     95.5       75 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.01245  0.1116  
 Residual               0.16744  0.4092  
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.167 

Conditional model:
                                            Estimate Std. Error z value
(Intercept)                                 -0.02766    0.07257  -0.381
cat_pre_wt_log_scale                        -3.64523    0.41435  -8.797
cat_pre_wt_log_scale_sq                     -3.01946    0.45628  -6.618
sl_mean_obs_log_scale                       -0.13590    0.05626  -2.416
prop_explore_logit_scale                     0.02370    0.04889   0.485
var_toxic_12_scale                          -0.01949    0.05115  -0.381
mean_toxic_conc_scale                       -0.25503    0.07147  -3.568
area_herb_log_scale                          0.76522    0.09726   7.868
beta_numeric_scale                          -0.11265    0.04837  -2.329
prop_explore_logit_scale:var_toxic_12_scale  0.18687    0.04488   4.164
                                            Pr(>|z|)    
(Intercept)                                 0.703071    
cat_pre_wt_log_scale                         < 2e-16 ***
cat_pre_wt_log_scale_sq                     3.65e-11 ***
sl_mean_obs_log_scale                       0.015710 *  
prop_explore_logit_scale                    0.627788    
var_toxic_12_scale                          0.703095    
mean_toxic_conc_scale                       0.000359 ***
area_herb_log_scale                         3.61e-15 ***
beta_numeric_scale                          0.019873 *  
prop_explore_logit_scale:var_toxic_12_scale 3.13e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.2 Variance in toxin (12 hour)

subm_var_toxic <- glmmTMB(
  var_toxic_12_scale ~ 
    sl_mean_obs_log_scale +
    beta_numeric_scale + 
    var_high + 
    (1|session_id),
  data = d2
); summary(subm_var_toxic)
 Family: gaussian  ( identity )
Formula:          
var_toxic_12_scale ~ sl_mean_obs_log_scale + beta_numeric_scale +  
    var_high + (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   178.2    193.0    -83.1    166.2       81 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 6.229e-11 7.892e-06
 Residual               3.956e-01 6.290e-01
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.396 

Conditional model:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)           -0.75075    0.09477  -7.922 2.34e-15 ***
sl_mean_obs_log_scale  0.21267    0.07751   2.744  0.00607 ** 
beta_numeric_scale    -0.15689    0.06811  -2.303  0.02126 *  
var_high               1.56021    0.13868  11.250  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.3 Mean step length

subm_sl <- glmmTMB(
  sl_mean_obs_log_scale ~ 
    prop_explore_logit_scale +  
    cat_pre_wt_log_scale + 
    cat_pre_wt_log_scale_sq + 
    on_toxic_logit_scale + 
    (1|session_id),
  data = d2,
); summary(subm_sl)
 Family: gaussian  ( identity )
Formula:          
sl_mean_obs_log_scale ~ prop_explore_logit_scale + cat_pre_wt_log_scale +  
    cat_pre_wt_log_scale_sq + on_toxic_logit_scale + (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   223.4    240.7   -104.7    209.4       80 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 4.610e-10 2.147e-05
 Residual               6.502e-01 8.064e-01
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.65 

Conditional model:
                         Estimate Std. Error z value Pr(>|z|)   
(Intercept)               0.06163    0.09646   0.639  0.52287   
prop_explore_logit_scale -0.26175    0.08987  -2.912  0.00359 **
cat_pre_wt_log_scale      1.37948    0.70121   1.967  0.04915 * 
cat_pre_wt_log_scale_sq   1.61095    0.75569   2.132  0.03303 * 
on_toxic_logit_scale      0.21435    0.09489   2.259  0.02389 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.4 Mean toxin ingested

subm_toxin_ingested <- glmmTMB(
  mean_toxic_conc_scale ~ 
    beta_numeric_scale + 
    on_toxic_logit_scale + 
    var_high + 
    (1|session_id),
  data = d2 
); summary(subm_toxin_ingested)
 Family: gaussian  ( identity )
Formula:          
mean_toxic_conc_scale ~ beta_numeric_scale + on_toxic_logit_scale +  
    var_high + (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   145.7    160.5    -66.8    133.7       81 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 6.461e-10 2.542e-05
 Residual               2.722e-01 5.217e-01
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.272 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)           0.18349    0.07861   2.334   0.0196 *  
beta_numeric_scale   -0.11683    0.05808  -2.012   0.0443 *  
on_toxic_logit_scale  0.28426    0.06393   4.447 8.72e-06 ***
var_high             -0.80646    0.11613  -6.945 3.79e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.5 Mean time on toxin

subm_on_toxic <- glmmTMB(
  on_toxic_logit_scale ~ 
    var_high + beta_numeric_scale * cat_pre_wt_log_scale + 
    (1|session_id),
  family = gaussian(),
  data = d2
); summary(subm_on_toxic)
 Family: gaussian  ( identity )
Formula:          
on_toxic_logit_scale ~ var_high + beta_numeric_scale * cat_pre_wt_log_scale +  
    (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   230.0    247.3   -108.0    216.0       80 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.01667  0.1291  
 Residual               0.68757  0.8292  
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.688 

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)  
(Intercept)                              0.22509    0.14140   1.592   0.1114  
var_high                                -0.41809    0.18032  -2.319   0.0204 *
beta_numeric_scale                      -0.12409    0.09766  -1.271   0.2039  
cat_pre_wt_log_scale                    -0.19306    0.13787  -1.400   0.1614  
beta_numeric_scale:cat_pre_wt_log_scale -0.23850    0.11732  -2.033   0.0421 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.6 Area consumption

subm_herb <- glmmTMB(
  area_herb_log_scale ~ 
    var_high * cat_pre_wt_log_scale +
    (1|session_id),
  data = d2,
); summary(subm_herb)
 Family: gaussian  ( identity )
Formula:          area_herb_log_scale ~ var_high * cat_pre_wt_log_scale + (1 |  
    session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   117.9    132.7    -53.0    105.9       81 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 3.501e-11 5.917e-06
 Residual               1.978e-01 4.448e-01
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.198 

Conditional model:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    0.04339    0.06963   0.623  0.53313    
var_high                       0.27072    0.10127   2.673  0.00751 ** 
cat_pre_wt_log_scale           0.55603    0.08378   6.637  3.2e-11 ***
var_high:cat_pre_wt_log_scale -0.26289    0.11558  -2.275  0.02293 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.7 Exploration

subm_exp <- glmmTMB(
  prop_explore_logit_scale ~ 
    var_high + beta_numeric_scale * cat_pre_wt_log_scale + 
    (1|session_id),
  family = gaussian(),
  data = d2
); summary(subm_exp)
 Family: gaussian  ( identity )
Formula:          
prop_explore_logit_scale ~ var_high + beta_numeric_scale * cat_pre_wt_log_scale +  
    (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   250.1    267.3   -118.0    236.1       80 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 3.781e-10 1.944e-05
 Residual               8.830e-01 9.397e-01
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.883 

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)  
(Intercept)                             -0.17231    0.14397  -1.197   0.2314  
var_high                                 0.15661    0.20355   0.769   0.4416  
beta_numeric_scale                      -0.14059    0.11040  -1.273   0.2029  
cat_pre_wt_log_scale                    -0.04358    0.12362  -0.353   0.7244  
beta_numeric_scale:cat_pre_wt_log_scale  0.25977    0.13242   1.962   0.0498 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.2 d-separation test

sem_fit <- psem(
    subm_rgr,
    subm_herb,
    subm_var_toxic,
    subm_on_toxic,
    subm_toxin_ingested,
    subm_exp,
    subm_sl, 
    var_toxic_12_scale %~~% on_toxic_logit_scale, # correlated errors
    data = d2
)
sem_summary <- suppressWarnings(suppressMessages(
  summary_psem(sem_fit, 
               no_standardize_x =  c(), 
               .progressBar = FALSE, 
               direction = "var_toxic_12_scale -> area_herb_log_scale")
))
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
sem_summary$Cstat
  Fisher.C df P.Value
1   42.469 42   0.451
knitr::kable(sem_summary$dTable)
Independ.Claim Test.Type DF Crit.Value P.Value
mean_toxic_conc_scale ~ cat_pre_wt_log_scale + … coef 87 -1.5712 0.1161
var_toxic_12_scale ~ cat_pre_wt_log_scale + … coef 87 -0.1892 0.8499
prop_explore_logit_scale ~ cat_pre_wt_log_scale_sq + … coef 87 0.8624 0.3884
mean_toxic_conc_scale ~ cat_pre_wt_log_scale_sq + … coef 87 1.8002 0.0718
area_herb_log_scale ~ cat_pre_wt_log_scale_sq + … coef 87 -0.7057 0.4804
on_toxic_logit_scale ~ cat_pre_wt_log_scale_sq + … coef 87 0.7944 0.4269
var_toxic_12_scale ~ cat_pre_wt_log_scale_sq + … coef 87 0.3081 0.7580
sl_mean_obs_log_scale ~ beta_numeric_scale + … coef 87 1.1020 0.2704
area_herb_log_scale ~ beta_numeric_scale + … coef 87 0.0200 0.9840
RGR_scale ~ var_high + … coef 87 0.6991 0.4845
sl_mean_obs_log_scale ~ var_high + … coef 87 -1.5883 0.1122
on_toxic_logit_scale ~ RGR_scale + … coef 87 -1.3664 0.1718
mean_toxic_conc_scale ~ sl_mean_obs_log_scale + … coef 87 0.7459 0.4557
area_herb_log_scale ~ sl_mean_obs_log_scale + … coef 87 0.8564 0.3918
mean_toxic_conc_scale ~ prop_explore_logit_scale + … coef 87 1.0591 0.2895
area_herb_log_scale ~ prop_explore_logit_scale + … coef 87 0.3920 0.6951
on_toxic_logit_scale ~ prop_explore_logit_scale + … coef 87 0.2361 0.8134
var_toxic_12_scale ~ prop_explore_logit_scale + … coef 87 1.1110 0.2666
area_herb_log_scale ~ mean_toxic_conc_scale + … coef 87 -0.2457 0.8059
var_toxic_12_scale ~ mean_toxic_conc_scale + … coef 87 -1.3086 0.1907
on_toxic_logit_scale ~ area_herb_log_scale + … coef 87 -0.7213 0.4707

3.3 Model summary

3.3.1 Model stats

sem_summary$ChiSq
NULL
sem_summary$AIC
       AIC     AICc  K  n
1 1264.813 1276.432 51 87
sem_summary$R2
                  Response   family     link method Marginal Conditional
1                RGR_scale gaussian identity   none     0.75        0.77
2      area_herb_log_scale gaussian identity   none     0.42          NA
3       var_toxic_12_scale gaussian identity   none     0.60          NA
4     on_toxic_logit_scale gaussian identity   none     0.19        0.21
5    mean_toxic_conc_scale gaussian identity   none     0.55          NA
6 prop_explore_logit_scale gaussian identity   none     0.06          NA
7    sl_mean_obs_log_scale gaussian identity   none     0.19          NA

3.3.2 Standardized coefficients

knitr::kable(
  cbind(sem_summary$coefficients,"arrow_size"=(abs(as.numeric(sem_summary$coefficients$Std.Estimate))*20))
)
Response Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate star arrow_size
RGR_scale cat_pre_wt_log_scale -3.6452 0.4144 87 -8.7974 0.0000 -3.5981 *** 71.962
RGR_scale cat_pre_wt_log_scale_sq -3.0195 0.4563 87 -6.6176 0.0000 -2.7708 *** 55.416
RGR_scale sl_mean_obs_log_scale -0.1359 0.0563 87 -2.4156 0.0157 -0.1448 * 2.896
RGR_scale prop_explore_logit_scale 0.0237 0.0489 87 0.4848 0.6278 0.0274 0.548
RGR_scale var_toxic_12_scale -0.0195 0.0511 87 -0.3811 0.7031 -0.0230 0.460
RGR_scale mean_toxic_conc_scale -0.2550 0.0715 87 -3.5683 0.0004 -0.2353 *** 4.706
RGR_scale area_herb_log_scale 0.7652 0.0973 87 7.8679 0.0000 0.5338 *** 10.676
RGR_scale beta_numeric_scale -0.1126 0.0484 87 -2.3287 0.0199 -0.1337 * 2.674
RGR_scale prop_explore_logit_scale:var_toxic_12_scale 0.1869 0.0449 87 4.1640 0.0000 0.2280 *** 4.560
area_herb_log_scale var_high 0.2707 0.1013 87 2.6733 0.0075 0.2315 ** 4.630
area_herb_log_scale cat_pre_wt_log_scale 0.5560 0.0838 87 6.6371 0.0000 0.7868 *** 15.736
area_herb_log_scale var_high:cat_pre_wt_log_scale -0.2629 0.1156 87 -2.2745 0.0229 -0.2795 * 5.590
var_toxic_12_scale sl_mean_obs_log_scale 0.2127 0.0775 87 2.7437 0.0061 0.1918 ** 3.836
var_toxic_12_scale beta_numeric_scale -0.1569 0.0681 87 -2.3033 0.0213 -0.1575 * 3.150
var_toxic_12_scale var_high 1.5602 0.1387 87 11.2503 0.0000 0.7875 *** 15.750
on_toxic_logit_scale var_high -0.4181 0.1803 87 -2.3186 0.0204 -0.2235 * 4.470
on_toxic_logit_scale beta_numeric_scale -0.1241 0.0977 87 -1.2706 0.2039 -0.1320 2.640
on_toxic_logit_scale cat_pre_wt_log_scale -0.1931 0.1379 87 -1.4003 0.1614 -0.1708 3.416
on_toxic_logit_scale beta_numeric_scale:cat_pre_wt_log_scale -0.2385 0.1173 87 -2.0329 0.0421 -0.2116 * 4.232
mean_toxic_conc_scale beta_numeric_scale -0.1168 0.0581 87 -2.0116 0.0443 -0.1503 * 3.006
mean_toxic_conc_scale on_toxic_logit_scale 0.2843 0.0639 87 4.4467 0.0000 0.3438 *** 6.876
mean_toxic_conc_scale var_high -0.8065 0.1161 87 -6.9447 0.0000 -0.5214 *** 10.428
prop_explore_logit_scale var_high 0.1566 0.2035 87 0.7694 0.4416 0.0808 1.616
prop_explore_logit_scale beta_numeric_scale -0.1406 0.1104 87 -1.2734 0.2029 -0.1443 2.886
prop_explore_logit_scale cat_pre_wt_log_scale -0.0436 0.1236 87 -0.3526 0.7244 -0.0372 0.744
prop_explore_logit_scale beta_numeric_scale:cat_pre_wt_log_scale 0.2598 0.1324 87 1.9617 0.0498 0.2224 * 4.448
sl_mean_obs_log_scale prop_explore_logit_scale -0.2618 0.0899 87 -2.9124 0.0036 -0.2840 ** 5.680
sl_mean_obs_log_scale cat_pre_wt_log_scale 1.3795 0.7012 87 1.9673 0.0491 1.2778 * 25.556
sl_mean_obs_log_scale cat_pre_wt_log_scale_sq 1.6109 0.7557 87 2.1318 0.0330 1.3873 * 27.746
sl_mean_obs_log_scale on_toxic_logit_scale 0.2144 0.0949 87 2.2589 0.0239 0.2244 * 4.488
~~var_toxic_12_scale ~~on_toxic_logit_scale 0.5323 - 87 5.7626 0.0000 0.5323 *** 10.646

3.3.3 Post hoc contrast

contrast_by_pre_wt(subm_herb, "var_high", "trt")
$emmeans
cat_pre_wt_log_scale = -2:
 var_high emmean    SE df lower.CL upper.CL
        0 -1.069 0.200 81   -1.467   -0.670
        1 -0.272 0.198 81   -0.666    0.122

cat_pre_wt_log_scale =  2:
 var_high emmean    SE df lower.CL upper.CL
        0  1.155 0.161 81    0.836    1.475
        1  0.900 0.150 81    0.603    1.198

Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast              estimate    SE df t.ratio p.value
 var_high1 - var_high0    0.797 0.281 81   2.830  0.0059

cat_pre_wt_log_scale =  2:
 contrast              estimate    SE df t.ratio p.value
 var_high1 - var_high0   -0.255 0.219 81  -1.162  0.2485
sd(subm_herb$frame$var_high)/sd(subm_herb$frame$area_herb_log_scale)
[1] 0.8549509
emmeans::emtrends(subm_exp, 
                 ~ beta_numeric_scale | cat_pre_wt_log_scale, 
                 var = "beta_numeric_scale", 
                 at = list(cat_pre_wt_log_scale = c(-2,2)))
cat_pre_wt_log_scale = -2:
 beta_numeric_scale beta_numeric_scale.trend    SE df lower.CL upper.CL
           7.88e-17                   -0.660 0.323 80    -1.30  -0.0173

cat_pre_wt_log_scale =  2:
 beta_numeric_scale beta_numeric_scale.trend    SE df lower.CL upper.CL
           7.88e-17                    0.379 0.246 80    -0.11   0.8677

Results are averaged over the levels of: var_high 
Confidence level used: 0.95 
sd(subm_exp$frame$beta_numeric_scale) * sd(subm_exp$frame$prop_explore_logit_scale)
[1] 0.974435
emmeans::emtrends(subm_on_toxic, 
                 ~ beta_numeric_scale | cat_pre_wt_log_scale, 
                 var = "beta_numeric_scale", 
                 at = list(cat_pre_wt_log_scale = c(-2,2)))
cat_pre_wt_log_scale = -2:
 beta_numeric_scale beta_numeric_scale.trend    SE df lower.CL upper.CL
           7.88e-17                    0.353 0.286 80   -0.217    0.923

cat_pre_wt_log_scale =  2:
 beta_numeric_scale beta_numeric_scale.trend    SE df lower.CL upper.CL
           7.88e-17                   -0.601 0.217 80   -1.034   -0.169

Results are averaged over the levels of: var_high 
Confidence level used: 0.95 
sd(subm_on_toxic$frame$beta_numeric_scale) * sd(subm_on_toxic$frame$on_toxic_logit_scale)
[1] 0.9403361

3.4 SEM plots

knitr::include_graphics(paste0(getwd(), "/graphs/manuscript1_figures/SEM_simplified.png"), rel_path = FALSE)

4 ISSF parameters

4.1 Selection (immigration) strength state 1

4.1.1 Model summary

m_select_s1_full <- glmmTMB(
  s1.less_toxic_end.est ~ 
    cat_pre_wt_log_scale * (beta + var_trt) + 
    (1|session_id),
  data = d3,
); summary(m_select_s1_full)
 Family: gaussian  ( identity )
Formula:          
s1.less_toxic_end.est ~ cat_pre_wt_log_scale * (beta + var_trt) +  
    (1 | session_id)
Data: d3

     AIC      BIC   logLik deviance df.resid 
   113.7    137.7    -46.9     93.7       71 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 3.344e-11 5.782e-06
 Residual               1.862e-01 4.315e-01
Number of obs: 81, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.186 

Conditional model:
                                     Estimate Std. Error z value Pr(>|z|)
(Intercept)                         -0.035699   0.098240  -0.363    0.716
cat_pre_wt_log_scale                 0.124186   0.114435   1.085    0.278
beta0                                0.095757   0.119051   0.804    0.421
beta5                               -0.036750   0.127318  -0.289    0.773
var_trtlow_var                       0.057820   0.100926   0.573    0.567
cat_pre_wt_log_scale:beta0           0.061099   0.137359   0.445    0.656
cat_pre_wt_log_scale:beta5           0.007268   0.146519   0.050    0.960
cat_pre_wt_log_scale:var_trtlow_var -0.078628   0.119191  -0.660    0.509
car::Anova(m_select_s1_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: s1.less_toxic_end.est
                              Chisq Df Pr(>Chisq)
(Intercept)                  0.1320  1     0.7163
cat_pre_wt_log_scale         1.1777  1     0.2778
beta                         1.2190  2     0.5436
var_trt                      0.3282  1     0.5667
cat_pre_wt_log_scale:beta    0.2245  2     0.8938
cat_pre_wt_log_scale:var_trt 0.4352  1     0.5095
m_select_s1 <- glmmTMB(
  s1.less_toxic_end.est ~ 
    cat_pre_wt_log_scale + beta + var_trt + 
    (1|session_id),
  data = d3,
); summary(m_select_s1)
 Family: gaussian  ( identity )
Formula:          
s1.less_toxic_end.est ~ cat_pre_wt_log_scale + beta + var_trt +  
    (1 | session_id)
Data: d3

     AIC      BIC   logLik deviance df.resid 
   108.3    125.1    -47.1     94.3       74 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 3.315e-11 5.758e-06
 Residual               1.875e-01 4.331e-01
Number of obs: 81, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.188 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)  
(Intercept)          -0.03455    0.09761  -0.354   0.7234  
cat_pre_wt_log_scale  0.10584    0.05788   1.829   0.0675 .
beta0                 0.11265    0.11575   0.973   0.3304  
beta5                -0.02728    0.11988  -0.228   0.8200  
var_trtlow_var        0.03783    0.09662   0.392   0.6954  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m_select_s1, type = "II")
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: s1.less_toxic_end.est
                      Chisq Df Pr(>Chisq)  
cat_pre_wt_log_scale 3.3438  1    0.06746 .
beta                 1.5799  2    0.45387  
var_trt              0.1533  1    0.69540  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.1.2 Post-hoc contrasts

contrast_by_pre_wt(m_select_s1, "cat_pre_wt_log_scale")
$emmeans
 cat_pre_wt_log_scale emmean    SE df lower.CL upper.CL
                   -2 -0.199 0.139 74 -0.47656   0.0788
                    2  0.225 0.112 74  0.00197   0.4471

Results are averaged over the levels of: beta, var_trt 
Confidence level used: 0.95 

$contrasts
 contrast                                         estimate    SE df t.ratio
 (cat_pre_wt_log_scale-2) - cat_pre_wt_log_scale2   -0.423 0.232 74  -1.829
 p.value
  0.0715

Results are averaged over the levels of: beta, var_trt 

4.1.3 Misc stats

r2(m_select_s1, tolerance = 10^-10)
Random effect variances not available. Returned R2 does not account for random effects.
# R2 for Mixed Models

  Conditional R2: NA
     Marginal R2: 0.057

4.1.4 Check model

check_mod(m_select_s1)

Warning: Maximum value of original data is not included in the
  replicated data.
  Model may not capture the variation of the data.

4.1.5 Plot

marginal_effects(m_select_s1, terms = c("cat_pre_wt_log_scale")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(yhat))) + 
  geom_hline(aes(yintercept = 1), color = "tomato", linetype = "dashed", linewidth = 1) + 
  geom_ribbon(aes(ymax = exp(upper), ymin = exp(lower)), alpha = 0.2) + 
  geom_line(linewidth = 2) + 
  geom_point(
    data = m_select_s1$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(s1.less_toxic_end.est)),
    size = 3, alpha = 0.8
  ) +  
  theme_bw(base_size = 15) + 
  theme(legend.position = "top") +
  scale_y_continuous(trans = "log10") + 
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  labs(x = "Pre-weight (g)", y = "Less toxic diet selection strength (odds ratio)")

4.2 Selection (immigration) strength state 2

4.2.1 Model summary

m_select_s2_full <- glmmTMB(
  s2.less_toxic_end.est ~ 
    cat_pre_wt_log_scale * (beta + var_trt) + 
    (1|session_id),
  data = d3 %>% 
      filter(s2.less_toxic_end.se < 10) # Filter out estimates that did not converge
); summary(m_select_s2_full)
 Family: gaussian  ( identity )
Formula:          
s2.less_toxic_end.est ~ cat_pre_wt_log_scale * (beta + var_trt) +  
    (1 | session_id)
Data: d3 %>% filter(s2.less_toxic_end.se < 10)

     AIC      BIC   logLik deviance df.resid 
   156.8    180.0    -68.4    136.8       65 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 9.311e-11 9.649e-06
 Residual               3.627e-01 6.023e-01
Number of obs: 75, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.363 

Conditional model:
                                     Estimate Std. Error z value Pr(>|z|)
(Intercept)                          0.145871   0.138823   1.051    0.293
cat_pre_wt_log_scale                 0.076065   0.163210   0.466    0.641
beta0                               -0.051651   0.168531  -0.306    0.759
beta5                               -0.064804   0.188723  -0.343    0.731
var_trtlow_var                      -0.006117   0.146957  -0.042    0.967
cat_pre_wt_log_scale:beta0           0.215500   0.194328   1.109    0.267
cat_pre_wt_log_scale:beta5          -0.241814   0.223029  -1.084    0.278
cat_pre_wt_log_scale:var_trtlow_var -0.178505   0.177311  -1.007    0.314
car::Anova(m_select_s2_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: s2.less_toxic_end.est
                              Chisq Df Pr(>Chisq)
(Intercept)                  1.1041  1     0.2934
cat_pre_wt_log_scale         0.2172  1     0.6412
beta                         0.1483  2     0.9285
var_trt                      0.0017  1     0.9668
cat_pre_wt_log_scale:beta    3.9324  2     0.1400
cat_pre_wt_log_scale:var_trt 1.0135  1     0.3141
# Refit
m_select_s2 <- glmmTMB(
  s2.less_toxic_end.est ~ 
    cat_pre_wt_log_scale + (var_trt +  beta) + 
    (1|session_id),
  data = d3 %>% filter(s2.less_toxic_end.se < 10) # Filter out estimates that did not converge
); summary(m_select_s2)
 Family: gaussian  ( identity )
Formula:          
s2.less_toxic_end.est ~ cat_pre_wt_log_scale + (var_trt + beta) +  
    (1 | session_id)
Data: d3 %>% filter(s2.less_toxic_end.se < 10)

     AIC      BIC   logLik deviance df.resid 
   154.8    171.0    -70.4    140.8       68 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 6.274e-11 7.921e-06
 Residual               3.827e-01 6.186e-01
Number of obs: 75, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.383 

Conditional model:
                      Estimate Std. Error z value Pr(>|z|)
(Intercept)           0.141025   0.141228   0.999    0.318
cat_pre_wt_log_scale  0.004266   0.086221   0.050    0.961
var_trtlow_var       -0.025702   0.144845  -0.177    0.859
beta0                 0.005924   0.166898   0.035    0.972
beta5                -0.119943   0.183479  -0.654    0.513
car::Anova(m_select_s2, type = "II")
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: s2.less_toxic_end.est
                      Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 0.0024  1     0.9605
var_trt              0.0315  1     0.8592
beta                 0.5632  2     0.7546

4.2.2 Post-hoc contrasts

emmeans::emmeans(m_select_s2, ~1)
 1       emmean     SE df lower.CL upper.CL
 overall 0.0913 0.0724 68  -0.0532    0.236

Results are averaged over the levels of: var_trt, beta 
Confidence level used: 0.95 

4.2.3 Misc stats

r2(m_select_s2, tolerance = 10^-10)
Random effect variances not available. Returned R2 does not account for random effects.
# R2 for Mixed Models

  Conditional R2: NA
     Marginal R2: 0.008

4.2.4 Check model

# Diagnostic plot looks better with m_select_s2 than m_select_s2_2
check_mod(m_select_s2)

Warning: Minimum value of original data is not included in the
  replicated data.
  Model may not capture the variation of the data.

4.2.5 Plot

marginal_effects(m_select_s2, terms = c("cat_pre_wt_log_scale")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(yhat))) + 
  geom_hline(aes(yintercept = 1), color = "tomato", linetype = "dashed", linewidth = 1) + 
  geom_ribbon(aes(ymax = exp(upper), ymin = exp(lower)), alpha = 0.2) + 
  geom_line(linewidth = 2) + 
  geom_point(
    data = m_select_s2$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(s2.less_toxic_end.est)),
    size = 3, alpha = 0.8
  ) +  
  theme_bw(base_size = 15) + 
  theme(legend.position = "top") +
  scale_y_continuous(trans = "log10") + 
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  labs(x = "Pre-weight (g)", y = "Less toxic diet selection strength (odds ratio)")

4.3 Arrestment state 1

4.3.1 Model summary

m_arrest_s1_full <- glmmTMB(
  log(k1) ~ 
    cat_pre_wt_log_scale * (beta + var_trt) + 
    (1|session_id),
  data = d3,
); summary(m_arrest_s1_full)
 Family: gaussian  ( identity )
Formula:          
log(k1) ~ cat_pre_wt_log_scale * (beta + var_trt) + (1 | session_id)
Data: d3

     AIC      BIC   logLik deviance df.resid 
   202.2    226.2    -91.1    182.2       71 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.000537 0.02317 
 Residual               0.554954 0.74495 
Number of obs: 81, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.555 

Conditional model:
                                    Estimate Std. Error z value Pr(>|z|)  
(Intercept)                         -0.08923    0.17091  -0.522   0.6016  
cat_pre_wt_log_scale                 0.04337    0.20406   0.212   0.8317  
beta0                                0.08266    0.20567   0.402   0.6877  
beta5                               -0.09266    0.21995  -0.421   0.6736  
var_trtlow_var                      -0.04815    0.17442  -0.276   0.7825  
cat_pre_wt_log_scale:beta0           0.10559    0.23857   0.443   0.6581  
cat_pre_wt_log_scale:beta5           0.63298    0.25306   2.501   0.0124 *
cat_pre_wt_log_scale:var_trtlow_var  0.09904    0.20584   0.481   0.6304  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m_arrest_s1_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: log(k1)
                              Chisq Df Pr(>Chisq)  
(Intercept)                  0.2726  1    0.60161  
cat_pre_wt_log_scale         0.0452  1    0.83170  
beta                         0.6378  2    0.72694  
var_trt                      0.0762  1    0.78251  
cat_pre_wt_log_scale:beta    6.7997  2    0.03338 *
cat_pre_wt_log_scale:var_trt 0.2315  1    0.63039  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m_arrest_s1 <- glmmTMB(
  log(k1) ~ 
    cat_pre_wt_log_scale * beta + var_trt + 
    (1|session_id),
  data = d3,
); summary(m_arrest_s1)
 Family: gaussian  ( identity )
Formula:          
log(k1) ~ cat_pre_wt_log_scale * beta + var_trt + (1 | session_id)
Data: d3

     AIC      BIC   logLik deviance df.resid 
   200.5    222.0    -91.2    182.5       72 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev.
 session_id (Intercept) 0.0001243 0.01115 
 Residual               0.5569498 0.74629 
Number of obs: 81, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.557 

Conditional model:
                           Estimate Std. Error z value Pr(>|z|)  
(Intercept)                -0.09552    0.17047  -0.560   0.5753  
cat_pre_wt_log_scale        0.09510    0.17338   0.548   0.5834  
beta0                       0.07431    0.20529   0.362   0.7174  
beta5                      -0.10004    0.21979  -0.455   0.6490  
var_trtlow_var             -0.02513    0.16800  -0.150   0.8811  
cat_pre_wt_log_scale:beta0  0.11796    0.23764   0.496   0.6196  
cat_pre_wt_log_scale:beta5  0.61691    0.25131   2.455   0.0141 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m_arrest_s1, type = "II")
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(k1)
                           Chisq Df Pr(>Chisq)   
cat_pre_wt_log_scale      6.8531  1   0.008849 **
beta                      0.3005  2   0.860488   
var_trt                   0.0224  1   0.881103   
cat_pre_wt_log_scale:beta 6.5451  2   0.037911 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m_arrest_s1, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: log(k1)
                           Chisq Df Pr(>Chisq)  
(Intercept)               0.3139  1    0.57527  
cat_pre_wt_log_scale      0.3008  1    0.58336  
beta                      0.6279  2    0.73057  
var_trt                   0.0224  1    0.88110  
cat_pre_wt_log_scale:beta 6.5451  2    0.03791 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.3.2 Post-hoc contrasts

contrast_by_pre_wt(m_arrest_s1, "beta")
$emmeans
cat_pre_wt_log_scale = -2:
 beta  emmean    SE df lower.CL upper.CL
 -5   -0.2983 0.409 72   -1.114    0.517
 0    -0.4599 0.432 72   -1.320    0.401
 5    -1.6321 0.492 72   -2.612   -0.652

cat_pre_wt_log_scale =  2:
 beta  emmean    SE df lower.CL upper.CL
 -5    0.0821 0.340 72   -0.597    0.761
 0     0.3923 0.370 72   -0.345    1.129
 5     1.2159 0.356 72    0.505    1.926

Results are averaged over the levels of: var_trt 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast         estimate    SE df t.ratio p.value
 (beta-5) - beta0    0.162 0.559 72   0.289  0.9550
 (beta-5) - beta5    1.334 0.613 72   2.178  0.0819
 beta0 - beta5       1.172 0.615 72   1.905  0.1446

cat_pre_wt_log_scale =  2:
 contrast         estimate    SE df t.ratio p.value
 (beta-5) - beta0   -0.310 0.473 72  -0.656  0.7894
 (beta-5) - beta5   -1.134 0.476 72  -2.381  0.0514
 beta0 - beta5      -0.824 0.483 72  -1.705  0.2104

Results are averaged over the levels of: var_trt 
Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 

4.3.3 Misc stats

r2(m_arrest_s1, tolerance = 10^-10)
# R2 for Mixed Models

  Conditional R2: 0.175
     Marginal R2: 0.175

4.3.4 Check model

check_mod(m_arrest_s1)

Warning: Maximum value of original data is not included in the
  replicated data.
  Model may not capture the variation of the data.

4.3.5 Plot

marginal_effects(m_arrest_s1, terms = c("cat_pre_wt_log_scale", "beta")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = (yhat), color = beta)) + 
  geom_ribbon(aes(ymax = (upper), ymin = lower, fill = beta, color = NULL), alpha = 0.2) + 
  geom_line(linewidth = 2) + 
  geom_point(
    data = m_arrest_s1$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(`log(k1)`)),
    size = 3, alpha = 0.5
  ) +  
  theme_bw(base_size = 15) + 
    geom_hline(aes(yintercept = 1), color = "black", linetype = "dashed", linewidth = 1) + 
  theme(legend.position = "top") +
  scale_y_continuous(trans = "log10") + 
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  scale_color_brewer(type = "qual", aesthetics = c("color","fill")) + 
  labs(x = "Pre-weight (g)", y = "Less toxic diet arresetment strength (ratio)", 
       color = expression(beta), fill = expression(beta))

4.4 Arrestment state 2

4.4.1 Model summary

m_arrest_s2_full <- glmmTMB(
  log(k2) ~ 
    cat_pre_wt_log_scale * (beta + var_trt) + 
    (1|session_id),
  data = d3,
); summary(m_arrest_s2_full)
 Family: gaussian  ( identity )
Formula:          
log(k2) ~ cat_pre_wt_log_scale * (beta + var_trt) + (1 | session_id)
Data: d3

     AIC      BIC   logLik deviance df.resid 
     154      178      -67      134       71 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.05334  0.2310  
 Residual               0.28154  0.5306  
Number of obs: 81, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.282 

Conditional model:
                                    Estimate Std. Error z value Pr(>|z|)  
(Intercept)                          0.20619    0.15939   1.294   0.1958  
cat_pre_wt_log_scale                -0.13519    0.18141  -0.745   0.4561  
beta0                               -0.09484    0.14743  -0.643   0.5201  
beta5                               -0.18572    0.15717  -1.182   0.2373  
var_trtlow_var                      -0.25005    0.12598  -1.985   0.0472 *
cat_pre_wt_log_scale:beta0           0.07897    0.16980   0.465   0.6419  
cat_pre_wt_log_scale:beta5           0.10713    0.18053   0.593   0.5529  
cat_pre_wt_log_scale:var_trtlow_var  0.09679    0.14916   0.649   0.5164  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m_arrest_s2_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: log(k2)
                              Chisq Df Pr(>Chisq)  
(Intercept)                  1.6735  1    0.19579  
cat_pre_wt_log_scale         0.5554  1    0.45612  
beta                         1.4088  2    0.49442  
var_trt                      3.9397  1    0.04716 *
cat_pre_wt_log_scale:beta    0.4036  2    0.81725  
cat_pre_wt_log_scale:var_trt 0.4211  1    0.51637  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m_arrest_s2 <- glmmTMB(
  log(k2) ~ 
    cat_pre_wt_log_scale + beta + var_trt + 
    (1|session_id),
  data = d3,
); summary(m_arrest_s2)
 Family: gaussian  ( identity )
Formula:          
log(k2) ~ cat_pre_wt_log_scale + beta + var_trt + (1 | session_id)
Data: d3

     AIC      BIC   logLik deviance df.resid 
   148.8    165.6    -67.4    134.8       74 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.04395  0.2096  
 Residual               0.28700  0.5357  
Number of obs: 81, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.287 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)  
(Intercept)           0.19258    0.15345   1.255   0.2095  
cat_pre_wt_log_scale -0.01354    0.13190  -0.103   0.9182  
beta0                -0.08219    0.14427  -0.570   0.5689  
beta5                -0.16897    0.14871  -1.136   0.2558  
var_trtlow_var       -0.23314    0.12060  -1.933   0.0532 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m_arrest_s2, type = "II")
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(k2)
                      Chisq Df Pr(>Chisq)  
cat_pre_wt_log_scale 0.0105  1    0.91823  
beta                 1.2923  2    0.52406  
var_trt              3.7367  1    0.05323 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.4.2 Post-hoc contrasts

emmeans::emmeans(m_arrest_s2, ~ "var_trt")
 var_trt  emmean    SE df lower.CL upper.CL
 high_var  0.105 0.132 74   -0.158    0.369
 low_var  -0.128 0.137 74   -0.400    0.145

Results are averaged over the levels of: beta 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

4.4.3 Misc stats

r2(m_arrest_s2, tolerance = 10^-10)
# R2 for Mixed Models

  Conditional R2: 0.176
     Marginal R2: 0.050

4.4.4 Check model

check_mod(m_arrest_s2)

Warning: Minimum value of original data is not included in the
  replicated data.
  Model may not capture the variation of the data.

4.4.5 Plot

marginal_effects(m_arrest_s2, terms = c("cat_pre_wt_log_scale", "beta")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = (yhat), color = beta)) + 
  geom_ribbon(aes(ymax = (upper), ymin = lower, fill = beta, color = NULL), alpha = 0.2) + 
  geom_line(linewidth = 2) + 
  geom_point(
    data = m_arrest_s2$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(`log(k2)`)),
    size = 3, alpha = 0.5
  ) +  
  theme_bw(base_size = 15) + 
    geom_hline(aes(yintercept = 1), color = "black", linetype = "dashed", linewidth = 1) + 
  theme(legend.position = "top") +
  scale_y_continuous(trans = "log10") + 
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  scale_color_brewer(type = "qual", aesthetics = c("color","fill")) + 
  labs(x = "Pre-weight (g)", y = "Less toxic diet arresetment strength (ratio)", 
       color = expression(beta), fill = expression(beta))

4.5 Toxic diet step length state 1

4.5.1 Model summary

m_sl_pred_s1_full <- glmmTMB(
  log(sl_mean_pred1) ~ 
    cat_pre_wt_log_scale * (beta + var_trt) + 
    (1|session_id),
  data = d3,
); summary(m_sl_pred_s1_full)
 Family: gaussian  ( identity )
Formula:          
log(sl_mean_pred1) ~ cat_pre_wt_log_scale * (beta + var_trt) +  
    (1 | session_id)
Data: d3

     AIC      BIC   logLik deviance df.resid 
   262.6    286.5   -121.3    242.6       71 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 8.455e-10 2.908e-05
 Residual               1.170e+00 1.082e+00
Number of obs: 81, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 1.17 

Conditional model:
                                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          3.98102    0.24627  16.165   <2e-16 ***
cat_pre_wt_log_scale                -0.03586    0.28686  -0.125    0.901    
beta0                                0.19843    0.29844   0.665    0.506    
beta5                                0.43478    0.31916   1.362    0.173    
var_trtlow_var                       0.34703    0.25300   1.372    0.170    
cat_pre_wt_log_scale:beta0           0.10587    0.34433   0.307    0.758    
cat_pre_wt_log_scale:beta5          -0.12070    0.36729  -0.329    0.742    
cat_pre_wt_log_scale:var_trtlow_var  0.22539    0.29879   0.754    0.451    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m_sl_pred_s1_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: log(sl_mean_pred1)
                                Chisq Df Pr(>Chisq)    
(Intercept)                  261.3177  1     <2e-16 ***
cat_pre_wt_log_scale           0.0156  1     0.9005    
beta                           1.8565  2     0.3952    
var_trt                        1.8814  1     0.1702    
cat_pre_wt_log_scale:beta      0.3619  2     0.8345    
cat_pre_wt_log_scale:var_trt   0.5690  1     0.4506    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m_sl_pred_s1 <- glmmTMB(
  log(sl_mean_pred1) ~ 
    cat_pre_wt_log_scale + beta + var_trt + 
    (1|session_id),
  data = d3,
); summary(m_sl_pred_s1)
 Family: gaussian  ( identity )
Formula:          log(sl_mean_pred1) ~ cat_pre_wt_log_scale + beta + var_trt +  
    (1 | session_id)
Data: d3

     AIC      BIC   logLik deviance df.resid 
   257.8    274.6   -121.9    243.8       74 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev.
 session_id (Intercept) 7.726e-10 2.78e-05
 Residual               1.188e+00 1.09e+00
Number of obs: 81, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 1.19 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)           3.95860    0.24562  16.117   <2e-16 ***
cat_pre_wt_log_scale  0.08645    0.14566   0.594   0.5528    
beta0                 0.20629    0.29126   0.708   0.4788    
beta5                 0.36034    0.30166   1.195   0.2323    
var_trtlow_var        0.41258    0.24315   1.697   0.0897 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m_sl_pred_s1, type = "II")
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: log(sl_mean_pred1)
                      Chisq Df Pr(>Chisq)  
cat_pre_wt_log_scale 0.3523  1    0.55282  
beta                 1.4497  2    0.48440  
var_trt              2.8793  1    0.08973 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.5.2 Post-hoc contrasts

emmeans::emmeans(m_sl_pred_s1, ~ "var_trt")
 var_trt  emmean    SE df lower.CL upper.CL
 high_var   4.17 0.172 74     3.83     4.51
 low_var    4.58 0.171 74     4.24     4.92

Results are averaged over the levels of: beta 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
emmeans::emmeans(m_sl_pred_s1, ~ 1)
 1       emmean    SE df lower.CL upper.CL
 overall   4.38 0.121 74     4.13     4.62

Results are averaged over the levels of: beta, var_trt 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

4.5.3 Misc stats

r2(m_sl_pred_s1, tolerance = 10^-10)
# R2 for Mixed Models

  Conditional R2: 0.053
     Marginal R2: 0.053

4.5.4 Check model

check_mod(m_sl_pred_s1)

4.5.5 Plot

marginal_effects(m_sl_pred_s1, terms = c("cat_pre_wt_log_scale")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = (yhat) / 1000 * 12)) + 
  geom_ribbon(aes(ymax = (upper)/ 1000 * 12, ymin = lower/ 1000 * 12), alpha = 0.2) + 
  geom_line(linewidth = 2) + 
  geom_point(
    data = m_sl_pred_s1$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), 
        y = exp(`log(sl_mean_pred1)`) / 1000 * 12),
    size = 3, alpha = 0.5
  ) +  
  theme_bw(base_size = 15) + 
  theme(legend.position = "top") +
  scale_y_continuous(trans = "log10", labels = fancy_scientific) + 
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  scale_color_brewer(type = "qual", aesthetics = c("color","fill")) + 
  labs(x = "Cat pre-weight (g)", y = "Selection free step length \non toxic diet (cm)")

4.6 Toxic diet step length state 2

4.6.1 Model summary

m_sl_pred_s2_full <- glmmTMB(
  log(sl_mean_pred2) ~ 
    cat_pre_wt_log_scale * (beta + var_trt) + 
    (1|session_id),
  data = d3,
); summary(m_sl_pred_s2_full)
 Family: gaussian  ( identity )
Formula:          
log(sl_mean_pred2) ~ cat_pre_wt_log_scale * (beta + var_trt) +  
    (1 | session_id)
Data: d3

     AIC      BIC   logLik deviance df.resid 
    71.5     95.5    -25.8     51.5       71 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 6.870e-12 2.621e-06
 Residual               1.106e-01 3.326e-01
Number of obs: 81, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.111 

Conditional model:
                                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          2.014673   0.075709  26.611   <2e-16 ***
cat_pre_wt_log_scale                 0.036250   0.088189   0.411   0.6810    
beta0                               -0.159917   0.091747  -1.743   0.0813 .  
beta5                               -0.083506   0.098117  -0.851   0.3947    
var_trtlow_var                      -0.068098   0.077779  -0.876   0.3813    
cat_pre_wt_log_scale:beta0           0.156273   0.105856   1.476   0.1399    
cat_pre_wt_log_scale:beta5          -0.131417   0.112915  -1.164   0.2445    
cat_pre_wt_log_scale:var_trtlow_var -0.008615   0.091854  -0.094   0.9253    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Somehow the Wald chisquare test is significant.
car::Anova(m_sl_pred_s2_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: log(sl_mean_pred2)
                                Chisq Df Pr(>Chisq)    
(Intercept)                  708.1363  1    < 2e-16 ***
cat_pre_wt_log_scale           0.1690  1    0.68103    
beta                           3.0394  2    0.21878    
var_trt                        0.7666  1    0.38128    
cat_pre_wt_log_scale:beta      6.2659  2    0.04359 *  
cat_pre_wt_log_scale:var_trt   0.0088  1    0.92528    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m_sl_pred_s2 <- glmmTMB(
  log(sl_mean_pred2) ~ 
    cat_pre_wt_log_scale * beta + var_trt + 
    (1|session_id),
  data = d3,
); summary(m_sl_pred_s2)
 Family: gaussian  ( identity )
Formula:          log(sl_mean_pred2) ~ cat_pre_wt_log_scale * beta + var_trt +  
    (1 | session_id)
Data: d3

     AIC      BIC   logLik deviance df.resid 
    69.5     91.1    -25.8     51.5       72 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 7.862e-12 2.804e-06
 Residual               1.106e-01 3.326e-01
Number of obs: 81, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.111 

Conditional model:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)                 2.01525    0.07547  26.704   <2e-16 ***
cat_pre_wt_log_scale        0.03168    0.07348   0.431   0.6664    
beta0                      -0.15920    0.09143  -1.741   0.0817 .  
beta5                      -0.08285    0.09787  -0.847   0.3973    
var_trtlow_var             -0.07011    0.07476  -0.938   0.3483    
cat_pre_wt_log_scale:beta0  0.15516    0.10519   1.475   0.1402    
cat_pre_wt_log_scale:beta5 -0.13003    0.11195  -1.162   0.2454    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m_sl_pred_s2, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: log(sl_mean_pred2)
                             Chisq Df Pr(>Chisq)    
(Intercept)               713.1060  1    < 2e-16 ***
cat_pre_wt_log_scale        0.1858  1    0.66641    
beta                        3.0330  2    0.21948    
var_trt                     0.8796  1    0.34830    
cat_pre_wt_log_scale:beta   6.4980  2    0.03881 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.6.2 Post-hoc contrasts

contrast_by_pre_wt(m_sl_pred_s2, "beta")
$emmeans
cat_pre_wt_log_scale = -2:
 beta emmean    SE df lower.CL upper.CL
 -5     1.92 0.174 72     1.57     2.26
 0      1.45 0.177 72     1.09     1.80
 5      2.09 0.209 72     1.68     2.51

cat_pre_wt_log_scale =  2:
 beta emmean    SE df lower.CL upper.CL
 -5     2.04 0.146 72     1.75     2.33
 0      2.19 0.150 72     1.90     2.49
 5      1.70 0.153 72     1.40     2.00

Results are averaged over the levels of: var_trt 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast         estimate    SE df t.ratio p.value
 (beta-5) - beta0    0.470 0.248 72   1.891  0.1486
 (beta-5) - beta5   -0.177 0.273 72  -0.650  0.7932
 beta0 - beta5      -0.647 0.274 72  -2.360  0.0540

cat_pre_wt_log_scale =  2:
 contrast         estimate    SE df t.ratio p.value
 (beta-5) - beta0   -0.151 0.209 72  -0.724  0.7503
 (beta-5) - beta5    0.343 0.212 72   1.616  0.2454
 beta0 - beta5       0.494 0.214 72   2.310  0.0608

Results are averaged over the levels of: var_trt 
Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans::emmeans(m_sl_pred_s2, ~ 1)
 1       emmean     SE df lower.CL upper.CL
 overall   1.91 0.0372 72     1.84     1.98

Results are averaged over the levels of: beta, var_trt 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

4.6.3 Misc stats

r2(m_sl_pred_s2, tolerance = 10^-10)
Random effect variances not available. Returned R2 does not account for random effects.
# R2 for Mixed Models

  Conditional R2: NA
     Marginal R2: 0.120

4.6.4 Check model

check_mod(m_sl_pred_s2)

4.6.5 Plot

marginal_effects(m_sl_pred_s2, terms = c("cat_pre_wt_log_scale", "beta")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), 
             y = (yhat) / 1000 * 12, 
             color = beta)) + 
  geom_ribbon(aes(ymax = (upper)/ 1000 * 12, ymin = lower/ 1000 * 12, 
                  fill = beta, color = NULL), 
              alpha = 0.2) + 
  geom_line(linewidth = 2) + 
  geom_point(
    data = m_sl_pred_s2$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), 
        y = exp(`log(sl_mean_pred2)`) / 1000 * 12),
    size = 3, alpha = 0.5
  ) +  
  theme_bw(base_size = 15) + 
  theme(legend.position = "top") +
  scale_y_continuous(trans = "log10") + 
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  scale_color_brewer(type = "qual", aesthetics = c("color","fill")) + 
  labs(x = "Pre-weight (g)", y = "Selection free step length \non toxic diet (cm)")

5 ISSF simulation

5.1 Observed means

5.1.1 Summary

# m <- brm(
#   on_toxic ~ 
#     beta * cat_pre_wt_log_scale + var_trt + 
#     (1|session_id),
#   family = Beta(),
#   data = d %>% 
#     filter(var_trt != "constant"), 
#   chains = 4, 
#   cores = 4, 
#   prior = c(
#     set_prior("normal(0,1.5)", class = "b"),
#     set_prior("normal(0,1)", class = "Intercept"),
#     set_prior("cauchy(0,0.5)", class = "sd"),
#     set_prior("gamma(0.01, 0.01)", class = "phi")
#   ), 
#   iter = 4000, 
#   control = list(adapt_delta = 0.95)
# )
# saveRDS(m, "invisible/fitted_models/on_toxic_brm.rds")
m <- readRDS("invisible/fitted_models/on_toxic_brm.rds")

m
 Family: beta 
  Links: mu = logit; phi = identity 
Formula: on_toxic ~ beta * cat_pre_wt_log_scale + var_trt + (1 | session_id) 
   Data: d %>% filter(var_trt != "constant") (Number of observations: 97) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Group-Level Effects: 
~session_id (Number of levels: 5) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.24      0.19     0.01     0.72 1.00     2231     3150

Population-Level Effects: 
                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                     -0.73      0.23    -1.19    -0.27 1.00     4467
beta0                          0.14      0.22    -0.29     0.56 1.00     7373
beta5                         -0.20      0.24    -0.66     0.28 1.00     6607
cat_pre_wt_log_scale           0.10      0.20    -0.28     0.49 1.00     4581
var_trtlow_var                 0.36      0.18     0.01     0.72 1.00     8461
beta0:cat_pre_wt_log_scale    -0.30      0.24    -0.76     0.18 1.00     5456
beta5:cat_pre_wt_log_scale    -0.72      0.29    -1.29    -0.16 1.00     5511
                           Tail_ESS
Intercept                      4001
beta0                          5574
beta5                          5527
cat_pre_wt_log_scale           5668
var_trtlow_var                 5568
beta0:cat_pre_wt_log_scale     5516
beta5:cat_pre_wt_log_scale     6012

Family Specific Parameters: 
    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi     4.72      0.66     3.51     6.09 1.00     6699     5717

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

5.1.2 Check model

brms:::plot.brmsfit(m)

brms::pp_check(m,ndraws = 100)

5.1.3 Plot

brms::conditional_effects(m,effects = c("cat_pre_wt_log_scale:beta"), re_formula = NA)[[1]] %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = estimate__, color = beta)) + 
  geom_ribbon(aes(ymax = upper__, ymin = lower__, color = NULL, fill = beta), alpha = 0.2) + 
  geom_line(size = 2) + 
  geom_point(
    data = m$data,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = on_toxic),
    size = 5, 
    alpha = 0.5
  ) + 
  theme_bw(base_size = 15) + 
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  scale_color_brewer(type = "qual", aesthetics = c("fill", "color")) + 
  labs(x = "Pre-weight (g)", y = "Proprtion of time on toxic diet", color = expression(beta), fill = expression(beta))

5.1.4 Observed effect sizes

# Compute observed effect size
epred_draws(m, expand.grid("cat_pre_wt_log_scale" = c(-2, 2), 
                       "beta" = c(-5, 5),
                       "var_trt" = c("low_var","high_var"),
                       "session_id" = "foo")) %>% 
  group_by(cat_pre_wt_log_scale, beta, draw) %>% 
  summarise(val = mean(val)) %>% 
  group_by(
    cat_pre_wt_log_scale, draw
  ) %>% 
  summarise(
    delta = diff(val)
  ) %>% 
  group_by(cat_pre_wt_log_scale) %>% 
  do(as.data.frame(t(summarise_vec(.$delta))))
# A tibble: 2 × 6
# Groups:   cat_pre_wt_log_scale [2]
  cat_pre_wt_log_scale   mean median     sd   lower   upper
                 <dbl>  <dbl>  <dbl>  <dbl>   <dbl>   <dbl>
1                   -2  0.280  0.284 0.151  -0.0251  0.563 
2                    2 -0.285 -0.284 0.0970 -0.479  -0.0991

5.2 Group means

out <- read_csv("simulation/move_rules_sim.csv")

out <- out %>% 
  mutate(
    scale = round(scale / 1000 * 12, 2),
    rss = as.factor(round(exp(rss), 2))
  )

out %>% 
  group_by(b, rss, scale, k) %>% 
  summarise(
    mean = mean(on_toxic), se = se(on_toxic), sd = sd(on_toxic), 
    lower = quantile(on_toxic, probs = 0.025), upper = quantile(on_toxic, probs = 0.975)
  ) %>% 
  mutate(
    z = (mean - 0.5) / se,
    P = pnorm(abs(z), lower.tail = FALSE) * 2
  ) %>% 
  knitr::kable()
b rss scale k mean se sd lower upper z P
-5 0.8 0.1 0.25 0.9039186 0.0018517 0.0370341 0.8151349 0.9671079 218.1331185 0.0000000
-5 0.8 0.1 1.00 0.5905470 0.0069139 0.1382788 0.3143606 0.8511489 13.0962854 0.0000000
-5 0.8 0.1 4.00 0.1331119 0.0052640 0.1052799 0.0029970 0.3729520 -69.6976278 0.0000000
-5 0.8 0.3 0.25 0.8059266 0.0017740 0.0354790 0.7432567 0.8761239 172.4549517 0.0000000
-5 0.8 0.3 1.00 0.5962987 0.0038047 0.0760946 0.4274476 0.7482517 25.3102497 0.0000000
-5 0.8 0.3 4.00 0.1527298 0.0033212 0.0664248 0.0459041 0.3127872 -104.5604328 0.0000000
-5 0.8 1.0 0.25 0.6634316 0.0011939 0.0238778 0.6153596 0.7083666 136.8900920 0.0000000
-5 0.8 1.0 1.00 0.5805495 0.0016245 0.0324908 0.5124126 0.6473776 49.5828564 0.0000000
-5 0.8 1.0 4.00 0.2898601 0.0024510 0.0490205 0.1928072 0.3846653 -85.7354275 0.0000000
-5 1 0.1 0.25 0.8788961 0.0023185 0.0463708 0.7831918 0.9590410 163.4202417 0.0000000
-5 1 0.1 1.00 0.4908342 0.0061322 0.1226433 0.2497253 0.7343656 -1.4947137 0.1349892
-5 1 0.1 4.00 0.0950150 0.0043984 0.0879685 0.0009990 0.3289960 -92.0750307 0.0000000
-5 1 0.3 0.25 0.7522952 0.0020909 0.0418177 0.6762737 0.8292208 120.6642601 0.0000000
-5 1 0.3 1.00 0.5039910 0.0037582 0.0751641 0.3644605 0.6543457 1.0619450 0.2882606
-5 1 0.3 4.00 0.1090909 0.0024322 0.0486447 0.0359141 0.2138611 -160.7200891 0.0000000
-5 1 1.0 0.25 0.5968731 0.0012855 0.0257093 0.5444306 0.6463786 75.3604269 0.0000000
-5 1 1.0 1.00 0.4997977 0.0017462 0.0349246 0.4315185 0.5704296 -0.1158482 0.9077728
-5 1 1.0 4.00 0.2229146 0.0019965 0.0399302 0.1458541 0.2987013 -138.7848297 0.0000000
-5 1.25 0.1 0.25 0.8291958 0.0029629 0.0592574 0.7012987 0.9380869 111.1070747 0.0000000
-5 1.25 0.1 1.00 0.4059590 0.0067969 0.1359388 0.1637862 0.6994755 -13.8357759 0.0000000
-5 1.25 0.1 4.00 0.0732118 0.0038774 0.0775487 0.0009990 0.2867383 -110.0696915 0.0000000
-5 1.25 0.3 0.25 0.6839660 0.0022279 0.0445578 0.6013237 0.7672827 82.5740974 0.0000000
-5 1.25 0.3 1.00 0.3934191 0.0035980 0.0719601 0.2627123 0.5456044 -29.6222231 0.0000000
-5 1.25 0.3 4.00 0.0847577 0.0019719 0.0394377 0.0229520 0.1689061 -210.5814170 0.0000000
-5 1.25 1.0 0.25 0.5268457 0.0013041 0.0260821 0.4764985 0.5764735 20.5854906 0.0000000
-5 1.25 1.0 1.00 0.4200100 0.0017184 0.0343685 0.3496503 0.4895854 -46.5485105 0.0000000
-5 1.25 1.0 4.00 0.1670554 0.0017809 0.0356186 0.0978272 0.2427822 -186.9499930 0.0000000
0 0.8 0.1 0.25 0.9154396 0.0020826 0.0416514 0.8289960 0.9800450 199.4841251 0.0000000
0 0.8 0.1 1.00 0.6122827 0.0083901 0.1678026 0.2497003 0.8991009 13.3827165 0.0000000
0 0.8 0.1 4.00 0.1549850 0.0080624 0.1612470 0.0009990 0.5434815 -42.7933403 0.0000000
0 0.8 0.3 0.25 0.8413237 0.0018827 0.0376546 0.7662338 0.9171329 181.2919850 0.0000000
0 0.8 0.3 1.00 0.6001873 0.0047834 0.0956689 0.4034965 0.7932817 20.9446032 0.0000000
0 0.8 0.3 4.00 0.1455420 0.0043286 0.0865713 0.0268981 0.3247502 -81.8881337 0.0000000
0 0.8 1.0 0.25 0.6926149 0.0014591 0.0291815 0.6363636 0.7472527 132.0114712 0.0000000
0 0.8 1.0 1.00 0.5832243 0.0022441 0.0448828 0.4984016 0.6674575 37.0851532 0.0000000
0 0.8 1.0 4.00 0.2501199 0.0030531 0.0610614 0.1368631 0.3816434 -81.8455512 0.0000000
0 1 0.1 0.25 0.8816234 0.0030176 0.0603529 0.7442308 0.9780220 126.4640036 0.0000000
0 1 0.1 1.00 0.5025924 0.0085678 0.1713558 0.1747752 0.8043706 0.3025761 0.7622129
0 1 0.1 4.00 0.1114461 0.0065858 0.1317168 0.0009990 0.4890360 -58.9983792 0.0000000
0 1 0.3 0.25 0.7911538 0.0024125 0.0482507 0.6923077 0.8841658 120.6837732 0.0000000
0 1 0.3 1.00 0.5023701 0.0050301 0.1006024 0.3096653 0.6893357 0.4711874 0.6375069
0 1 0.3 4.00 0.1084366 0.0030472 0.0609447 0.0169081 0.2468282 -128.4979711 0.0000000
0 1 1.0 0.25 0.6282917 0.0016147 0.0322930 0.5664086 0.6933067 79.4547584 0.0000000
0 1 1.0 1.00 0.5003621 0.0021468 0.0429356 0.4135864 0.5884116 0.1686888 0.8660415
0 1 1.0 4.00 0.1918981 0.0024640 0.0492807 0.0948801 0.2918082 -125.0394879 0.0000000
0 1.25 0.1 0.25 0.8398027 0.0036314 0.0726281 0.6822927 0.9620629 93.5733645 0.0000000
0 1.25 0.1 1.00 0.4036014 0.0086244 0.1724887 0.1177073 0.7423576 -11.1773828 0.0000000
0 1.25 0.1 4.00 0.0906968 0.0053304 0.1066081 0.0009990 0.4089411 -76.7865320 0.0000000
0 1.25 0.3 0.25 0.7233292 0.0028702 0.0574040 0.6082667 0.8372128 77.8096042 0.0000000
0 1.25 0.3 1.00 0.4058666 0.0048732 0.0974646 0.2227023 0.5925325 -19.3164192 0.0000000
0 1.25 0.3 4.00 0.0789535 0.0023950 0.0479005 0.0129620 0.2027972 -175.8004437 0.0000000
0 1.25 1.0 0.25 0.5525574 0.0017095 0.0341894 0.4844156 0.6173826 30.7449156 0.0000000
0 1.25 1.0 1.00 0.4180020 0.0021056 0.0421124 0.3356643 0.4956543 -38.9424151 0.0000000
0 1.25 1.0 4.00 0.1420055 0.0020180 0.0403607 0.0629121 0.2238262 -177.3975346 0.0000000
5 0.8 0.1 0.25 0.9130170 0.0044098 0.0881967 0.6540210 0.9990010 93.6581027 0.0000000
5 0.8 0.1 1.00 0.5673002 0.0143028 0.2860567 0.0029471 1.0000000 4.7053742 0.0000025
5 0.8 0.1 4.00 0.2545105 0.0151281 0.3025619 0.0009990 1.0000000 -16.2273900 0.0000000
5 0.8 0.3 0.25 0.9005519 0.0025916 0.0518327 0.7900100 0.9820180 154.5557517 0.0000000
5 0.8 0.3 1.00 0.5884416 0.0087466 0.1749316 0.2246004 0.8841159 10.1115560 0.0000000
5 0.8 0.3 4.00 0.1503721 0.0071233 0.1424664 0.0009990 0.4257493 -49.0821650 0.0000000
5 0.8 1.0 0.25 0.7920380 0.0024577 0.0491536 0.7022727 0.8971029 118.8266235 0.0000000
5 0.8 1.0 1.00 0.5995904 0.0040889 0.0817777 0.4374875 0.7643856 24.3563692 0.0000000
5 0.8 1.0 4.00 0.1726299 0.0044209 0.0884190 0.0239011 0.3627373 -74.0497284 0.0000000
5 1 0.1 0.25 0.8940559 0.0048582 0.0971640 0.6282717 1.0000000 81.1114992 0.0000000
5 1 0.1 1.00 0.4758841 0.0137388 0.2747767 0.0069181 0.9960040 -1.7553080 0.0792066
5 1 0.1 4.00 0.1965135 0.0129869 0.2597375 0.0009990 0.9639361 -23.3687092 0.0000000
5 1 0.3 0.25 0.8604745 0.0033270 0.0665396 0.7150849 0.9680320 108.3489209 0.0000000
5 1 0.3 1.00 0.4881194 0.0095941 0.1918820 0.1195804 0.8563686 -1.2383254 0.2155954
5 1 0.3 4.00 0.1007742 0.0057848 0.1156956 0.0009990 0.4299451 -69.0131201 0.0000000
5 1 1.0 0.25 0.7348826 0.0029112 0.0582234 0.6383616 0.8551698 80.6832908 0.0000000
5 1 1.0 1.00 0.4977522 0.0041455 0.0829107 0.3036963 0.6643606 -0.5422101 0.5876738
5 1 1.0 4.00 0.1246628 0.0034988 0.0699754 0.0259241 0.2897103 -107.2769047 0.0000000
5 1.25 0.1 0.25 0.8516084 0.0062154 0.1243090 0.5332667 0.9980020 56.5700727 0.0000000
5 1.25 0.1 1.00 0.4237737 0.0129779 0.2595572 0.0108891 0.9451049 -5.8735634 0.0000000
5 1.25 0.1 4.00 0.1613786 0.0128102 0.2562041 0.0009990 0.9990260 -26.4337184 0.0000000
5 1.25 0.3 0.25 0.8178372 0.0041144 0.0822890 0.6378871 0.9641359 77.2490361 0.0000000
5 1.25 0.3 1.00 0.4056618 0.0086395 0.1727904 0.1228521 0.7922328 -10.9193729 0.0000000
5 1.25 0.3 4.00 0.0920504 0.0051176 0.1023529 0.0009990 0.3383866 -79.7142890 0.0000000
5 1.25 1.0 0.25 0.6606144 0.0032006 0.0640114 0.5593906 0.8052448 50.1830249 0.0000000
5 1.25 1.0 1.00 0.4090260 0.0040690 0.0813807 0.2657343 0.5804446 -22.3576409 0.0000000
5 1.25 1.0 4.00 0.0917183 0.0027186 0.0543716 0.0119630 0.2268981 -150.1820280 0.0000000

5.3 Some fun

# Precompute some summarize grid
z <- out %>% 
  filter(b != 0) %>% 
  filter(scale == 1) %>% 
  group_by(b, rss, k) %>%
  dplyr::select(on_toxic) %>% 
  mutate(
    rss = as.numeric(as.character(rss)),
    id = seq_along(on_toxic)
  )

# Observed draws
a <- epred_draws(m, expand.grid("cat_pre_wt_log_scale" = c(-2, 2), 
                           "beta" = c(-5, 5),
                           "var_trt" = c("low_var","high_var"),
                           "session_id" = "foo")) %>% 
  group_by(cat_pre_wt_log_scale, beta, draw) %>% 
  summarise(val = mean(val)) %>% 
  group_by(
    cat_pre_wt_log_scale, draw
  ) %>% 
  summarise(
    delta = diff(val)
  ) %>% 
  group_split(cat_pre_wt_log_scale, .keep = FALSE) %>% 
  lapply(
    function(x){
      x$delta
    }
  )

names(a) <- c("small", "large")



# Some simple wrapper for hypothesis testing later in the model scenarios
compute_delta <- function(hypothesis, res, ref_list){
  temp <- hypothesis %>% 
    left_join(res, by = c("k", "rss","b"), relationship = "many-to-many") %>% 
    group_by(cat_size, id) %>% 
    summarise(
      delta = diff(on_toxic)
    )
  out1 <- temp %>% 
    group_by(cat_size) %>% 
    do(as.data.frame(t(summarise_vec(.$delta)))) %>% 
    arrange(rev(cat_size))
  
  out2 <- temp %>% 
    group_split(cat_pre_wt_log_scale, .keep = FALSE) %>% 
    suppressMessages() %>% 
    suppressWarnings() %>% 
    lapply(function(x){
      x$delta
    })
  
  names(out2) <- c("large", "small")
  
  r1 <- bayestestR::overlap(ref_list$small, out2$small)
  r2 <- bayestestR::overlap(ref_list$large, out2$large)
  return(
    list(
      "effects" = out1,
      "overlap_coef" = sigfig(c("small" = r1, "large" = r2))
    )
  )
}

5.4 Model scenarios

# No arrestment, no immigration
data.frame("k" = c(1, 1, 1, 1), 
           "rss" = c(1, 1, 1, 1),
           "cat_size" = c("s", "s", "l", "l"),
           "b" = c(-5, 5, -5, 5)) %>% 
  compute_delta(z, a) 
$effects
# A tibble: 2 × 6
# Groups:   cat_size [2]
  cat_size     mean  median     sd  lower upper
  <chr>       <dbl>   <dbl>  <dbl>  <dbl> <dbl>
1 s        -0.00205 0.00300 0.0891 -0.176 0.183
2 l        -0.00205 0.00300 0.0891 -0.176 0.183

$overlap_coef
 small  large 
"0.23" "0.13" 
# Size dependent arrestment, no immigration
data.frame("k" = c(0.25, 0.25, 4, 4), 
           "rss" = c(1,1, 1, 1),
           "cat_size" = c("s", "s", "l", "l"),
           "b" = c(-5, 5, -5, 5)) %>% 
  compute_delta(z, a)
$effects
# A tibble: 2 × 6
# Groups:   cat_size [2]
  cat_size    mean median     sd   lower  upper
  <chr>      <dbl>  <dbl>  <dbl>   <dbl>  <dbl>
1 s         0.138   0.138 0.0639  0.0189 0.261 
2 l        -0.0983 -0.107 0.0803 -0.217  0.0783

$overlap_coef
 small  large 
"0.44" "0.28" 
# No arrestment, size dependent immigration
data.frame("k" = c(1, 1, 1, 1), 
           "rss" = c(0.8,0.8, 1.25, 1.25),
           "cat_size" = c("s", "s", "l", "l"),
           "b" = c(-5, 5, -5, 5)) %>% 
  compute_delta(z, a)
$effects
# A tibble: 2 × 6
# Groups:   cat_size [2]
  cat_size    mean   median     sd  lower upper
  <chr>      <dbl>    <dbl>  <dbl>  <dbl> <dbl>
1 s         0.0190  0.0210  0.0880 -0.166 0.192
2 l        -0.0110 -0.00849 0.0865 -0.185 0.155

$overlap_coef
 small  large 
"0.26" "0.15" 
# Size dependent arrestment, size dependent immigration
data.frame("k" = c(0.25, 0.25, 4, 4), 
           "rss" = c(0.8,0.8, 1.25, 1.25),
           "cat_size" = c("s", "s", "l", "l"),
           "b" = c(-5, 5, -5, 5)) %>% 
  compute_delta(z, a)
$effects
# A tibble: 2 × 6
# Groups:   cat_size [2]
  cat_size    mean  median     sd   lower  upper
  <chr>      <dbl>   <dbl>  <dbl>   <dbl>  <dbl>
1 s         0.129   0.125  0.0568  0.0339 0.247 
2 l        -0.0753 -0.0804 0.0651 -0.188  0.0850

$overlap_coef
 small  large 
"0.38" "0.18" 
# Size dependent arrestment in clustered treatment, no immigration
data.frame("k" = c(1, 0.25, 1, 4), 
           "rss" = c(1,1, 1, 1),
           "cat_size" = c("s", "s", "l", "l"),
           "b" = c(-5, 5, -5, 5)) %>% 
  compute_delta(z, a)
$effects
# A tibble: 2 × 6
# Groups:   cat_size [2]
  cat_size   mean median     sd  lower  upper
  <chr>     <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1 s         0.235  0.235 0.0693  0.104  0.374
2 l        -0.375 -0.380 0.0790 -0.500 -0.194

$overlap_coef
 small  large 
"0.62" "0.58" 
# Size dependent arrestment in clustered treatment, size dependent immigration
data.frame("k" = c(1, 0.25, 1, 4), 
           "rss" = c(0.8,0.8, 1.25, 1.25),
           "cat_size" = c("s", "s", "l", "l"),
           "b" = c(-5, 5, -5, 5)) %>% 
  compute_delta(z, a)
$effects
# A tibble: 2 × 6
# Groups:   cat_size [2]
  cat_size   mean median     sd  lower  upper
  <chr>     <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1 s         0.211  0.208 0.0592  0.110  0.331
2 l        -0.328 -0.332 0.0655 -0.445 -0.166

$overlap_coef
 small  large 
"0.53" "0.68" 

5.5 Plots

5.5.1 Simplified plot

out <- read_csv("simulation/move_rules_sim.csv")

out <- out %>% 
  mutate(
    scale = round(scale / 1000 * 12, 2),
    rss = as.factor(round(exp(rss), 2))
  )

out %>% 
  filter(b != 0) %>% 
  filter(scale == 1) %>% 
  group_by(b, rss, k) %>%
  dplyr::select(on_toxic) %>% 
  mutate(
    id = seq_along(on_toxic)
  ) %>% 
  arrange(b, rss, k, id) %>% 
  group_by(rss, k, id) %>% 
  summarise(
    delta = diff(on_toxic)
  ) %>% 
  ggplot(aes(x = as.factor(k), y = delta, color = factor(rss))) +
  geom_hline(aes(yintercept = 0), color = "grey", 
             size = 1, linetype = "dashed") + 
  geom_point(position = position_jitterdodge(jitter.width = 0.2, 
                                             jitter.height = 0, 
                                             dodge.width = 0.5), 
             alpha = 0.2, 
             size = 2) + 
  geom_point(stat = "summary", 
             position = position_dodge(width = 0.5),
             shape = "—",
             size = 4,
             color = "black",
             aes(group = factor(rss))) + 
  theme_bw(base_size = 15) +
  theme(legend.position = "top") + 
  guides(colour = guide_legend(
    override.aes = list(alpha = 1, size = 4), 
    title.position = "top"
  )) + 
  labs(x = "Less toxic diet arresetment strength (ratio)", 
       y = expression(atop(Change~"in"~proprtion~time~on, paste("more toxic diet",~(beta[5]-beta[-5])))),
       color = "Less toxic diet relative immigration strength (odds ratio)") + 
  scale_color_brewer(type = "seq") -> g;g

5.5.2 Full result

out <- read_csv("simulation/move_rules_sim.csv")

out <- out %>% 
  mutate(
    scale = round(scale / 1000 * 12, 2),
    rss = as.factor(round(exp(rss), 2))
  )

k_lab <- unique(out$k)
names(k_lab) = sprintf("k = %s", k_lab)
scale_lab <- unique(out$scale)
names(scale_lab) = sprintf("scale = %s", scale_lab)
rss_lab <- unique(out$rss)
names(rss_lab) = sprintf("rss = %s", rss_lab)


out %>% 
  ggplot(aes(x = as.factor(k), y = on_toxic, color = factor(b))) + 
  geom_point(position = position_jitterdodge(jitter.width = 0.2, 
                                             jitter.height = 0, 
                                             dodge.width = 0.5), 
             alpha = 0.2,
             size = 1) + 
  geom_point(stat = "summary", 
             position = position_dodge(width = 0.5),
             shape = "—",
             size = 3,
             color = "black",
             aes(group = factor(b))) + 
  theme_bw(base_size = 15) +
  guides(colour = guide_legend(
    override.aes = list(alpha = 1, size = 3), 
  )) + 
  facet_grid(scale ~ rss, labeller = labeller(
    rss = reverse_names(rss_lab), scale = reverse_names(scale_lab)
  )) + 
  guides(colour = guide_legend(override.aes = list(alpha = 1))) + 
  labs(x = "Less toxic diet arresetment strength (ratio)", 
       y = "Proportion time on toxic diet", 
       color = expression(beta)) + 
  scale_color_brewer(type = "qual")->g;g

6 Session Info

sessionInfo()
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] spat1f_0.0.1       herbivar_0.2.1     forcats_0.5.1      stringr_1.5.0     
 [5] dplyr_1.1.2        purrr_1.0.1        readr_2.1.2        tidyr_1.3.0       
 [9] tibble_3.2.1       tidyverse_1.3.1    extraDistr_1.9.1   imager_0.42.13    
[13] magrittr_2.0.3     doSNOW_1.0.20      snow_0.4-4         iterators_1.0.14  
[17] foreach_1.5.2      glmmTMB_1.1.7      tidybayes_3.0.2    performance_0.10.4
[21] sjPlot_2.8.15      ggpubr_0.6.0       ggplot2_3.5.1      piecewiseSEM_2.3.0
[25] survival_3.3-1    

loaded via a namespace (and not attached):
  [1] fs_1.5.2                 matrixStats_1.2.0        bitops_1.0-7            
  [4] sf_1.0-8                 devtools_2.4.5           lubridate_1.8.0         
  [7] httr_1.4.3               RColorBrewer_1.1-3       insight_0.19.3          
 [10] doParallel_1.0.17        numDeriv_2016.8-1.1      profvis_0.3.7           
 [13] tools_4.3.1              backports_1.4.1          DT_0.23                 
 [16] sjlabelled_1.2.0         utf8_1.2.2               R6_2.5.1                
 [19] mgcv_1.8-40              ggdist_3.1.1             urlchecker_1.0.1        
 [22] withr_2.5.0              sp_1.5-0                 readbitmap_0.1.5        
 [25] gridExtra_2.3            Brobdingnag_1.2-9        prettyunits_1.1.1       
 [28] cli_3.6.1                shinyjs_2.1.0            sandwich_3.0-2          
 [31] labeling_0.4.2           sass_0.4.1               mvtnorm_1.1-3           
 [34] spatstat.data_3.0-0      brms_2.19.0              proxy_0.4-27            
 [37] StanHeaders_2.26.13      DHARMa_0.4.6             MuMIn_1.47.5            
 [40] sessioninfo_1.2.2        readxl_1.4.0             rstudioapi_0.13         
 [43] circular_0.5-0           visNetwork_2.1.0         generics_0.1.2          
 [46] gtools_3.9.2.2           crosstalk_1.2.0          vroom_1.5.7             
 [49] car_3.1-2                distributional_0.3.0     inline_0.3.19           
 [52] loo_2.5.1                Matrix_1.5-4.1           fansi_1.0.3             
 [55] abind_1.4-5              lifecycle_1.0.3          multcomp_1.4-25         
 [58] yaml_2.3.5               carData_3.0-5            grid_4.3.1              
 [61] qgam_1.3.4               promises_1.2.0.1         crayon_1.5.1            
 [64] miniUI_0.1.1.1           lattice_0.21-8           haven_2.5.0             
 [67] cowplot_1.1.1            pillar_1.9.0             knitr_1.43              
 [70] boot_1.3-28.1            estimability_1.3         shinystan_2.6.0         
 [73] codetools_0.2-19         imagerExtra_2.0.0        glue_1.6.2              
 [76] RcppArmadillo_0.11.2.0.0 V8_4.2.0                 remotes_2.4.2           
 [79] vctrs_0.6.3              png_0.1-7                Rdpack_2.3.1            
 [82] cellranger_1.1.0         gtable_0.3.0             assertthat_0.2.1        
 [85] cachem_1.0.6             xfun_0.39                rbibutils_2.2.8         
 [88] mime_0.12                coda_0.19-4              shinythemes_1.2.0       
 [91] units_0.8-0              DiagrammeR_1.0.9         ellipsis_0.3.2          
 [94] fitdistrplus_1.1-11      TH.data_1.1-1            gap_1.2.3-6             
 [97] nlme_3.1-162             CircStats_0.2-6          xts_0.13.1              
[100] usethis_2.1.6            bit64_4.0.5              threejs_0.3.3           
[103] rstan_2.26.13            rprojroot_2.0.3          tensorA_0.36.2          
[106] bslib_0.3.1              TMB_1.9.4                KernSmooth_2.23-20      
[109] colorspace_2.0-3         DBI_1.1.3                tidyselect_1.2.0        
[112] processx_3.6.1           emmeans_1.7.5            curl_4.3.2              
[115] bit_4.0.4                compiler_4.3.1           amt_0.2.1.0             
[118] rvest_1.0.2              arrayhelpers_1.1-0       see_0.7.1               
[121] xml2_1.3.3               desc_1.4.1               colourpicker_1.1.1      
[124] bayestestR_0.13.1        posterior_1.2.2          dygraphs_1.1.1.6        
[127] checkmate_2.1.0          scales_1.3.0             classInt_0.4-7          
[130] callr_3.7.0              tiff_0.1-11              digest_0.6.29           
[133] fftwtools_0.9-11         spatstat.utils_3.0-1     minqa_1.2.4             
[136] rmarkdown_2.24           base64enc_0.1-3          htmltools_0.5.2         
[139] moveHMM_1.9              pkgconfig_2.0.3          jpeg_0.1-9              
[142] lme4_1.1-29              highr_0.9                dbplyr_2.2.0            
[145] fastmap_1.1.0            rlang_1.1.1              htmlwidgets_1.5.4       
[148] shiny_1.7.1              farver_2.1.0             jquerylib_0.1.4         
[151] svUnit_1.0.6             zoo_1.8-12               jsonlite_1.8.8          
[154] bayesplot_1.10.0         geosphere_1.5-18         munsell_0.5.0           
[157] Rcpp_1.0.8.3             stringi_1.7.12           MASS_7.3-60             
[160] plyr_1.8.7               pkgbuild_1.3.1           parallel_4.3.1          
[163] sjmisc_2.8.9             deldir_1.0-6             ggeffects_1.2.3         
[166] splines_4.3.1            hms_1.1.1                sjstats_0.18.2          
[169] ps_1.7.1                 igraph_1.3.2             spatstat.geom_3.0-3     
[172] markdown_1.1             ggsignif_0.6.3           reshape2_1.4.4          
[175] stats4_4.3.1             pkgload_1.3.2            rstantools_2.2.0        
[178] reprex_2.0.1             evaluate_0.15            RcppParallel_5.1.5      
[181] modelr_0.1.8             bmp_0.3                  nloptr_2.0.3            
[184] tzdb_0.3.0               httpuv_1.6.5             RgoogleMaps_1.4.5.3     
[187] polyclip_1.10-0          broom_1.0.5              gap.datasets_0.0.5      
[190] xtable_1.8-4             e1071_1.7-11             rstatix_0.7.2           
[193] later_1.3.0              class_7.3-22             memoise_2.0.1           
[196] ggmap_3.0.2              bridgesampling_1.1-2